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ABSTRACT 


A fast  minimum  mean-square  error  technique  for 
restoring  images  degraded  by  blur  is  presented  in  this 
dissertation . 


The  phenomenon  of  linear  shift-invariant  degradation 
in  an  incoherent  optical  system  can  be  described  by  means 
of  a convolution  integral.  Since  digital  processing  of 
pictorial  data  requires  discretization  of  this  integral  by 
means  of  a quadrature  technique,  a theoretical  study  of  a 
broad  class  of  quadrature  formulae  is  first  presented. 


The  discrete  image  degradation  phenomenon  is  modelled 
by  two  distinct  vector  opace  formulations:  dark  background 
objects  correspond  to  a model  possessing  an  overdetermined 
blur  matrix;  objects  with  unknown  background,  however, 
result  in  a system  that  is  underdetermined . It  is  shown 
that  these  models  become  equivalent  if  the  background  of 
the  object  is  artificially  set  to  zero  by  processing  the 
observed  image.  This  fact  results  in  introduction  of  a 
fast  restoration  technique  in  the  absence  of  noise. 


The  noisy  restoration  problem  is  resolved  by  employing 
Wiener  estimation.  It  is  shown  that  with  proper  arrangement 
of  the  observed  image  data,  the  covariance  matrix  of  the 
object  becomes  a circulant  matrix.  Hence,  the  Fourier 
domain  properties  of  circulants  gives  rise  to  a 


computationally  efficient  Wiener  restorat ^ f'.i  technique.  A 
certain  approximation  imposed  oi.  this  technique  results  in 
a suboptimal,  but  faster,  restoration  filter.  It  is  shown 
that  the  computational  saving  gained  by  this  approximation 
is  significant,  while  the  increase  in  the  error  variance  is 


quite  small. 
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1.  INTRODUCTION 


The  concept  of  digital  restoration  in  this  text  is 
interpreted  as  the  reconstruction  of  an  image  to  an 
original  object  by  the  removal  of  degradation  phenomenon 
known  a priori,  possibly  in  the  presence  of  noise  [1-1] . 
Thus,  restoration  techniques  require  some  form  of  knowledge 
concerning  the  degradation  phenomenon,  and  this  knowledge 
either  comes  from  a deterministic  assumption  about  the 
phenomenon  or  statistical  models.  Degradation  systems  are 
often  quite  complex  in  general.  However,  in  many  cases  of 
practical  importance,  the  degrading  systems  can  be 
presented  by  a linear  smoothing  operation  followed  by  the 
addition  of  noise  known  only  in  a statistical  sense  [1-2]. 


Early  image  restoration  techniques,  which  were 
introduced  by  the  pioneers  of  this  field  [1-3]  [1-4],  did 
not  prove  very  successful  because  the  proposed  techniques 
did  not  acknowledge  the  existence  of  noise  in  the  imaging 
systems.  To  be  more  specific,  the  failure  of  the  early 
restoration  methods  was  caused  because  of  the  amplification 
of  a high  frequency  noise  component  by  the  inversion  of 
some  degrciding  process.  Tsujiuchi  [1-5]  , Harris  [1-6]  , 
HcGlamery  [1-7] , and  Mueller  [1-8]  have  all  tried  to 
overcome  this  hindrance  by  modifying  the  inversion 
technique.  For  example,  Harris  [1-6]  has  replaced  the 
inverse  filter  by  zero  for  the  range  of  spatial  frequencies 
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tor  which  the  noise  power  exceeds  the  power  of  the  signsl. 
Lster,  researchers  in  this  field  realized  that  a mote 
successful  recovery  of  the  image  could  be  achieved  with  a 
„re  realistic  approach  which  utilized  the  characteristics 
of  the  noise  in  the  imaging  systems.  Consequently,  an 
optimum  .inear  shift-invariant  filter  was  introduced  by 
Helstrom  11-91  which,  when  applied  on  a noisy  degraded 
image,  gives  an  estimate  of  the  ideal  image  with  the  least 
mean-square  error.  This  filter  is  indeed  the  same  as  the 
classical  Wiener  filter.  Slepian  (1-101  has  also  solved 
the  same  problem  when  the  smoothing  function  itself  in 
stochastic.  It  has  later  be.in  illustrated  that  the  minimum 
mean-square  error  technique  gives  better  reconstructed 
images  than  the  simple  inversion  method  [1-111. 

The  Wiener  filter  technique  unf or  tuna tely  has  the 
limitations  of  large  storage  requirements  along  with 

inefficient  computational  methods.  Pratt  11-121  has 

Mionfzr  filtering  computation 
introduced  generalized  Wiener  tiireui.  j 

techniques  which,  by  utilizing  transform  properties  of 

imaging  systems,  have  improved  the  computational 

efficiency.  Furthermore,  Pratt  illustrated  that  a specific 

computational  procedure  could  result  in  a significant 

reduction  of  the  computational  load,  with  only  a small 

increase  in  estimation  error.  It  has  also  been  shown  that 

lower-triangular  transformations  can  give  rise  to  an 

efficient  suboptimal  Wiener  filter  [1-13]  . 
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Constrained  restoration  is  an  alternative  solution  for 
the  problem  of  noisy  image  reconstruction.  Here  the 
mean-square  error  is  not  necessarily  minimized,  but  high 
frequency  noise  oscillations  are  dampened  by  observing  the 
constraints  governing  the  image  forming  systems.  Hunt 
[1-14]  has  employed  special  properties  of  linear  systems  to 
introduce  fast  constrained  image  estimation  techniques.  A 
constrained  restoration  technique  introduced  by  Mascarenhas 
[1-15]  utilizes  linear  equality  and  inequality  constraints. 
Linear  inequality  constraints  involve  solution  of  a 
quadratic  programming  problem,  and  require  extensive 
computing  when  images  of  reasonable  size  are  to  be 
processed.  A specific  case  of  inequality  constrained  image 
restoration  is  when  positiveness  of  image  intensities  is 
utilized  for  better  image  reconstruction  purposes.  A 
survey  of  positive  image  restoration  techniques  is  given  in 
reference  [1-16]  . 


Blind  deconvolution,  in  which  the  point-spread 
function  is  assumed  unknown  [1-17] , has  attracted  some 
attention.  Ekstrom  [1-18]  has  suggested  means  of 
estimating  the  unknown  point-spread  function  by  processing 
the  blurred  and  noisy  observation.  Cole  [1-19]  has 
introduced  the  homomorphic  filter  which  is  the  geometrical 
mean  between  the  Wiener  filter  and  the  inverse  filter. 
Similar  methods  have  been  applied  to  the  problem  of 
restoring  old  acoustic  recordings  as  well  as  reconstruction 
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o£  blurry  Imaqes  11-19).  (For  a comparison  o£  different 
restoration  techniques,  the  reader  is  referred  to  reference 

[1-20] ) . 

The  fundamental  purpose  of  this  dissertation  is  to 
explore  restoration  techniques  which  ace  computationally 
fast  and  efficient,  but  can  be  applied  to  the  restoration 
of  large  sire  images.  This  work  starts  with  a brief,  but 
sufficiently  broad,  discussion  of  the  discretization 
methods  of  continuous  linear  shift-invariant  degradation 
systems.  Next,  the  problem  of  a fast  pseudoinverse 

technique  which  can  be  applied  to  a general  noise-free 
image  Is  resolved:  this  part  can  be  considered  to  be  an 
application  and  extension  of  the  Fourier  transform 

properties  of  clcculant  matrices  H-211 , (1-22).  In  the 

presence  of  noise,  a fast  Wiener  filtering  technique  is 
developed  which  overcomes  many  of  the  computational 
obstacles  of  the  conventional  Wiener  filter.  Finally,  the 
problem  of  fast  constrained  filtering  of  degraded  Images  is 

considered. 
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2.  THE  RESTORATION  PROBLEM  IN  ITS  CONTINUOUS  FORM 


2.  1 The  Image 


Every  visible  object  can  be  characterized  by  its 
radiant  energy  distribution  which  commonly  is  described  by 
a two- Umensional  function  of  spatial  variables,  f(x,y). 
The  two-dimensional  radiant  energy  distribution  wh_ch 
enters  the  human  eye  is  often  referred  to  as  the  image. 
Even  in  the  case  of  an  observer  with  perfect  vision,  the 
image  itself  might  be  a distorted  replica  of  the  object 
function.  This,  inevitably,  will  cause  an  imperfect 
comprehension  of  the  scene  by  the  brain.  A degraded  image 
can  result  from  many  different  phenomena.  A turbulent 
atmosphere,  for  example,  deforms  the  phase  function  while 
the  light  travels  tnrough  the  air.  When  collected  by  this 
optical  system  aperture,  a blurred  image  results.  A 
photographic  camera  with  a poor  lens  produces  a low  quality 
picture.  Also  a misfocused  lens  and  movements  of  the 
objects  in  a scene  both  generate  errors  in  recording  the 
scene,  and  thus  the  image  often  differs  from  the  original 
object  in  one  way  or  another.  The  following  section 
studies  the  mathematical  model  of  the  image-forming  and 
image-degrading  processes.  Since  t!iis  dissertation  deals 
exclusively  with  linear  phenomena,  only  degradations  which 
can  be  represented  by  a linear  system  have  been  considered. 
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2.  2 Convolution  as  a Model  for  Linear  Systems 


A continuous  space  invariant  ( 
system  is  interpreted  to  be  a convoluti 
a fixed  function  h,  commonly  referred 
response  or  the  point  spread  functi 
which  is  the  input  to  the  system, 
expression  for  convolution  is 
oo 


g (x)  = 


f (s) n (x-s) ds 


stationary)  linear 
on  of  two  functions: 
to  as  the  impulse 
on,  and  a function  f 
The  mathematical 


(2-1) 


where  g is  the  output  of  the  system.  In  a more  compact 
notation,  eq.  (2-1)  can  be  written  as 


g(x)=f(x)Wh(x) 


(2-2) 


The  output  of  a linear  system  is  equal  to  the  impulse 
response,  h,  when  the  input  to  the  system  is  an  impulse. 


In  two  dimensions  the  signals  f,  h,  and  g are 
functions  of  two  variables  as  modelled  by  the  relation 
Oo  oo 

g(x,y)=  / / f (r ,s) h (x-r  ,y-s) drds  (2-3) 


x,y)=  j j f (r ,s) h(x-r ,y-s) 


-CO  -oo 


The  two  dimensional  convolution  integral 
model  for  a spatially  invariant  degradation 
incoherent  illumination  [2-1],  [2-2]. 


is  the  proper 
occurring  under 


A linear  system  is  fully  defined  by  its  impulse 


response , 


And,  this  is  also  true  for  a linear  degradation 


where  the  impulse  response  of  a particular  process 
completely  characterizes  the  degradation  phenomenon. 


A diffraction-limited  rectangular  optical  system  is 
cliaracter  ized  by  a separable  impulse  response  of  the  form 
[2-3J 


2 2 
i sin  (r)»  (sin(  s)  i 
h(r,s)=  


(2-4) 


Blurring  due  to  atmospheric  turbulence  has  been  modelled  by 
a linear  operator  of  the  form  [2-4] 


2 2 

h(r ,s)=exp[- (r  +s  ) ) 


(2-5) 


Motion  blur  has  a one  dimensional  point  spread  function 
defined  as  [2-5] 


h's)=l 


if  -1/2^  s <1/2 


(2-6) 


otherwise 


Although  the  above  is  not  an  exhaustive  list  of  sources  of 
linear  image  degradation,  the  list  provides  a general 
intuition  for  some  of  the  problems  faced  in  image 
restoration.  The  next  section  describes  the  use  of  the 
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Fourier  transform  of  linear  systems  for  removal  of 
spatially  invariant  degradations. 


2.  3 Inverse  Filtering 


To  av^id  notational  complexity  and  unnecessary 
formulations,  only  one-dimensional  degradation  is  discussed 
at  this  point.  This  approach  will  be  followed  in  most 
remaining  sections  of  this  dissertation,  and  except  for  the 
examples  explicitly  discussed,  the  extension  of  a 
one-dimensional  problem  to  higher  dimensions  is  assumed  to 
be  straight  forward.  Considering  the  model  for  image 
degradation  formulated  by  eq.  (2-1),  the  restoration  task 
is  phrased  as  follows:  assuming  the  observation,  g,  and  the 
point  spread  function,  h,  are  both  given,  attempt  to 
recover  or  estimate  the  image,  f. 


It  has  been  shown  that  Fourier  techniques  can  play  an 
important  role  in  attempting  to  obtain  the  object  from  the 
observation  g(x)  through  the  inversion  of  h(.)  [2-6], 
[2-7].  For  a function  f(x)  the  Fourier  transform,  F(w),  of 
f(x)  is  defined  as 


F(w)  = 


L 


oo 

f (x) exp[-iwx] dx 


(2-7) 


The  above  integral  does  not,  however,  exist  for  every 
function  f(x)  [2-8],  but  the  existence  of  this  integral  for 
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the  class  of  functions  encountered  in  this  dissertation  is 
certain . 

An  interesting  utilization  of  the  i’ourier  transform  is 
its  application  to  linear  systems.  Let  G(w),  F(w),  and 
H(w)  denote  the  Fourier  transforms  of  g(x),  f(x),  and  h(x), 
respectively.  Then  from  eguation  (2-1)  it  is  easily  shown 
that 


G(w)=F(w)  H(w) 


(2-8) 


Observing  the  above  equality,  the  concept  of  inverse 
filtering  becomes  clear.  Inverse  filtering  simply  consist 
of  dividing  both  sides  or  eq.  (2-8)  by  the  blur  transfer 
function  H(w).  If  H(w)  does  not  vanish  at  any  point,  the 
object,  f(x),  can  be  completely  recovered  by  the  operation 


-1 

f(x)='^ 


G(w) 

HK) 


^ -1 

where  the  operation  ^ denotes  the 
transform 


(2-9) 

inverse  Fourier 


1 

f (x)  = ( ) 

2U 


oo 


F(w) expliwx] dw 


(2-10) 


Equation  (2-10)  is  the  inverse  of  eq.  {2-1),  and  the 
notation  is  the  same  in  both  equations. 
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Although  the  problem  of  noisy  observations  will  be 
studied  in  later  chapters,  a brief  coMent  on  eq.  (2-9) 
seems  essential  at  this  point.  The  impulse  response  H(w) 
almost  always  dec-eases  rapidly  with  growth  of  w.  On  the 
other  hand,  any  small  amount  of  noise  or  uncertainty  in  the 
Observation  has  a relatively  flat  spectral  distribution. 
This  means  that  the  inverse  filtering  technique  enhances 
high  frequency  noise  so  strongly  that  even  for  a small 
amount  of  noise  the  technique  could  not  be  applied,  and 
thus  must  be  modified. 
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3.  DISCRETIZATION  OF  THE  CONTINUOUS  MODEL 

To  treat  degraded  images  by  means  of  digital 
computers,  the  continuous  model  of  eq.  (2-1)  must  be 
discretized.  The  discretization  process  will,  naturally, 
replace  the  integral  by  a discrete  summation  which  takes 
advantage  of  values  of  the  signals  only  at  discrete  points. 

3.  1 Practical  Considerations 

In  most  practical  situations,  the  physical  sample 
image  g(x)  of  eg.  (2-1)  is  not  available  over  the  entire 
real  line,  and  axso  the  whole  infinite  extent  of  the  object 

is  not  usually  of  particular  interest  to  the  image 

processer.  In  fact,  in  practice,  only  finite  size  of 
objects  are  to  be  restored  by  processing  of  finite  size 
observations.  The  proceeding  argument  implies  that,  in 
reality,  the  limits  on  the  definite  integral  of  eq.  (2-1) 
are  not  infinite.  Another  important  feature  in  the  model 

of  eg.  (2-1)  is  that  the  degrading  function  h usually 

vanishes  beyond  some  point,  and  consequently  the  region  for 
which  h is  nonzero  has  a finite  length.  In  theory,  of 
course,  most  point  spread  functions  have  infinite  length, 
but  invariably,  h decreases  rapidly  for  large  values  of  x 
(the  examples  in  Sec.  (2-2),  for  instance,  have  this 
property).  Considering  this  characteristic,  a point  spread 
function  can  be  truncated  to  some  lengh  L without  severe 


16 


«dellinq  error  providing  that  the  length  L is  seiected 
wieely.  Figure  (3-1)  shows  a truncated  two-di.ensional 
Gaussian  point  spread  function. 


In  view  of  the  above  argument,  eg.  (2-1) 


is  modified 


to  the  form 


.v=x+-2 


g(x)=  I f(s)h(x-s)ds 


(3-1) 


' L 
u=x-^ 


where  [u,v]  is  the  region  of  integration. 


3.  2 Quadrature  formulae 

By  definition,  a quadrature  formula  (q.f.) 
approximation  to  a definite  integral,  the  approximation  is 
a linear  combination  of  values  of  the  integrand  and  its 
derivatives  at  certain  points  of  the  interval  of 
integration  called  the  nodes  of  the  q.f.  1 3-1).  When  the 
derivatives  of  the  integrand  are  unknown,  then  the  general 
form  of  a q.f.  is  expressed  as 


f (x)dx 


n 


)+Rf 


(3-2) 


where  c.  and  x.  are  the  coefficients  and  nodes  of  the  o.f.. 


respectively,  and  u<x,  <x^ <*„«•  "he  term  Rf  is  a 

functional  which,  for  any  given  function  f(*),  equals  the 
difference  between  the  exact  value  of  the  integral  and  its 


approximation,  and  Rf  may  vanish  for  some  class  of 
functions.  If  the  nodes  of  the  q.f.  are  selected  in 
advance,  the  only  available  parameters  to  be  treated  are 
the  coefficients  c.  . Examples  of  this  type  (fixed  node) 
are  the  method  called  pulse  approximation  (equal 
coefficient),  Newton-Cotes,  and  the  best  q.f.  in  the  sense 
of  Sard.  If  the  nodes  of  a q.f.  are  free,  the  best 
location  of  the  nodes,  in  a certain  sense,  can  be 
determined,  and  the  q.f.  is  called  optimal.  Examples  of 
the  optimal  type  are  Gauss-Legendre  and  the  optimal  q.f.  in 
the  sense  of  Sard.  In  image  restoration  the  nodes  x^  are 
usually  preassigned,  therefore,  only  fixed  node  quadrature 
formulae  are  discussed  here. 


3.  3 Spline  Functions  and  Sard's  Best  Quadrature  Formulae 


Given  a set  of  real  numbers 


X-=U  <x,<x_<....  <x  <v=x  .1 
0 12  n n+1 


(3-3) 


a spline  S(x)  of  degree  m with  the  nodes  x ,x  ,...,x  is  a 
function  defined  on  the  real  line  so  that  in  each  interval 


(X.  ,x.  , ),  for  i=0,l,...,n,  S(x)  is  represented  by  a 
polynomial  of  degree  m or  less,  and  the  function  and  its 
derivatives  of  order  m-1  or  less  are  continuous  on  [u,v]  , 
thus  S(x)  is  m-1  times  continuously  differentiable  [3-2]  . 


To  represent  splines,  truncated  power  functions  can  be 
employed  to  construct  a set  of  basis  functions  for  the 
spline  space.  A truncated  power  function  is  defined  as 


(3-4) 


A spline  function  of  deqree  m and  number  of  nodes  n, 

S (x) , has  a unique  representation  (3-2]  of  the  form 
m,  n 


m n 

S a x' + (mnY^c  (x-x.  ) 

m,  n i ^ ^ ^ + 

i = 0 i=l 


(3-5) 


where  a and  c are  unknown  coefficients  to  be  determined, 
i i 

TO  develop  Sard's  best  q.f.,  let  K(x)  be  a monospline 
of  degree  m [3-3]  with  n preassigned  nodes.  By  definition, 
a monospline  of  degree  m is  a spline  of  deqree  m-1  plus  a 
polynomial  of  deqree  m;  thus,  K(x)  can  be  formulated  as 


K(x)=— ‘S  (X) 

ml  m-l,n 


(3-6) 


It  is  known  that  an  arbitrary  monospline  can  give  rise  to  a 
q.f.  [3-4].  To  achieve  this,  set 


(i)  (i) 

K (u)=K  (v)=0 


(3-7) 


for  i=0,l,... ,m-l,  and  using  K(x)  as  the  kernel  [3-4],  then 
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n 


V 

^f(x)dx=^^c-  f(Xj)+Rf 


(3-8) 


where 


Rf=  I f (x)K(x)dx 


(3-9) 


Note 
less . 
norm) 


q.  f . 


that  Rf  vanishes  if  f is  a polynomial  of  degree  m 1 or 
If  K(x)  has  the  least  square  deviation  (minimum 
among  all  kernels  of  the  form  (3-6),  then  the 
(3-8)  is  called  best  in  the  sense  of  Sard.  Thus, 

V 


K 


2 

[K(x)  1 


dx=minimum 


(3-lfl) 


Schoenberg  [3-3],  [3-51  has  shown  that  there  exists 
unique  monospline  H(x)  of  degree  2m 


2m 

X 

H(x)= + S (x) 

(2m) ! 2m-l,n 


(3-11) 


with  nodes  Xj^,X2»..* 
three  conditions; 


,Xyj,  that  satisfies  the  following 


H(Xj )=0 

i — 

(m+i) 

H (u)=0 

(3-12a) 


(3-12b) 
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(3-120 


H (v)=0  fii'— 1 


In  terms  of  H(x),  the  kernel  K(x)  of  Sard  s best  q.f. 


given  by 


(m) 

K(x)=H  (X) 


(3-13) 


and  the  minimum  norm  of  K(x)  is  obtained  as 


■II-/' 

•'ll 


[ K(X)  f dx=  (-1)”^!  H(x)dx 


(3-14) 


By  normalizing  (u,v]  to  [-1,11  and  applying  condition 


(3-12b),  H(x)  simplifies  to 


H(x)  = 


n , „ .2m- 1 

(X4-1)  ^ ^ 


^ i=l 


(3-15) 


^ (2m-l)  ! 


conditions  (3-12a)  and  (3-120  produce  mtn  independent 
eouations,  whose  solution  gives  the  coefficients  a,  and  C;  . 
An  upper  bound  can  be  derived  for  the  error  ter.n  Rf  using 

eo.  (3-9).  Thus, 


r (m)  j (ni)| 

Rf=  I K(x)  f (x)dx<||K  jl  f 


(3-16) 


(3-17) 


1 

(m)|! 

Rf<(-1)"^/  H(x)dx  1 

f 11 

The  q.f.  of  form  eq.  (3-8)  with  coefficients  c 
obtained  from  eq . (3-12)  has  some  interesting  properties. 
By  varying  m from  1 to  n , eq.  (3-8)  presents  a large  family 
of  quadrature  formulae.  The  case  m=l,  if  x are  placed 
uniformly,  is  sometimes  called  tne  pulse  approximation 
method.  An  upper  bound  for  the  error  term  Rf  can  be 
derived  when  m<n.  This  property  is  an  important  one 

because  it  makes  study  of  the  error  possible  even  in  the 
simple  case  of  pulse  approximation.  When  m=n  the  technique 
is  called  Newton-Cotes  method.  Newton-Cotes  q.f.  results 
in  zero  error  for  the  class  of  polynomials  of  degree  n-1  or 
less,  but  no  explicit  error  term  is  given  if  the  integrand 
does  not  belong  to  this  class.  The  following  example  is 
designed  to  aid  the  reader  in  better  understanding  of  the 
best  q.f.  in  the  sense  of  Sard. 

Let  m=l  and  assume  Xj  are  uniformly  placed  on  [-1,1]. 
Figure  (3-2)  illustrates  the  location  of  the  nodes  for  the 
case  of  n=5 . When  the  nodes  are  placed  uniformly  on 
[-1,1],  the  location  of  the  nodes  is  obtained  from 


(3-18) 
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Here  the  expression  for  H(x)  is  given  by 


2 n 


H(x) 


(x+1) 

a c.  (x-x 

O 0 i i- 


(3-19) 


Equation  (3-12b)  then  reduces  to 


c =2 
f i 


(3-20) 


and  eq.  (3-12a)  becomes 


E 


(Xi  +1) 


Ci  (X  i“X;  ) -ao=- 


(3-21) 


For  i*l,2,...,n.  From  equations  (3-20)  and  (3-21),  and 

Ci  are  obtained  as 


(3--22a) 


Ci=  — 


(3-22b) 


for  i=l,2,...,n.  Substituting  (3-22b)  in  (3-8),  the  latter 
equation  becomes 


f (X) dx 


)+Rf 


(3-23) 


TO  estimate  Rf,  the  norm  of  K(x)  must  be  established. 


Thorefore , 


/ 


1 

H (x) dx=- 


(x+l)3 


n 


1-1 


2n 


(3-24) 


1-1  i=l 


4 1 (4n-l)  2 

3 n 3n  3n 


Thus  the  norm  of  K(x)  can  be  obtained  as 


H(x)dx=- 


3n 


(3-25) 


Using  the  inequality  of  form  eq.  (3-16) , the  upper  bound  of 
Rf  is  then  established  as 


Rf< 


I I 


3n 


And,  as  one  would  expect 


Lim  Rf=0 
n-«c>o 


(3-26) 


(3-27) 


In  the  process  of  determining  a q.f.,  one  is  faced 

with  the  task  of  selecting  parameter  m.  Equation  (3-9) 

(m^ 

requires  that,  for  a given  m,  f must  exist,  and  this 
condition  is  met  by  most  functions  dealt  with  in  this 
dissertation  for  any  value  of  m.  The  stability  of  the 
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discrete  system  as  well  as  the  error  term  Rf  are  affected 
by  the  choice  of  m.  Reference  [3-61  studies  the  behavior 
of  Rf  for  different  classes  of  functions  regarding  specific 
values  of  m.  In  view  of  the  conclusion  at  (3-ul,  m=l,  2, 
and  3 are  considered  to  be  good  choices  for  the  problems 
discussed  in  this  dissertation. 


3.  4 The  Overdetermined  Model 

A discrete  image  degradation  system  which  is  of  full 
column  rank  is  characterized  as  overdetermined.  In 
practice  this  situation  arises  from  either  a dark 
background  in  the  object  scene,  or.  over-sampling  of  the 
observation  [3-7].  Assuming  that  the  object  exists  only  on 
some  given  interval  [a,b],  the  observation  at  a given  point 
X is  formulated  as 


•'a 


f (s)  h(x^-s)ds 


where  x-  is  in  [a rb+ — ) and  h(x) 

^ 2 2 

function.  Thus, 


(3-28) 


is  a space  limited 


h (X. -s) =0 


if  X.  -s  > — 


(3-29) 


Considering  the  above  condition,  eg.  (3-28)  becomes 
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x+L/2 

g(x.  )=f  (s)h(x  ,-s)ds 


(^-30) 


X.-  L/2 
1 


Using  a q.f.  on  the  integral  of  ea.  (3-30)  and  taking  the 
nodes  at  a discrete  set  of  points,  the  above  model  is 
discretized  as 


g(x,)-^c.  f(s,)h(x.-s.: 


(3-31) 


If  the  nodes  s are  placed  uniformly  (a  valid  assumption  in 
J 

regard  to  image  restoration  problems) , s ^ can  be  assumed  to 
coincide  with  the  integers  on  the  real  line  without  loss  of 
generality.  Thus 


g (X.  ) . f(j)h(x.-j) 


(3-32) 


Equation  (3-32)  formulates  a general  discrete  image 
degradation  process.  If  the  observation  g(x)  is  sampled 
uniformly  at  points  x^=i,  then 


•xL-1 


q(i)  = > c.  f(j)h(i-j) 


(3-33) 


. .,L-1 


Employing  vector  space  notation,  eq.  (3-33)  can  be 
presented  in  a more  compact  form.  To  construct  the  vector 
space  model,  assume  the  object  is  defined  with  M samples  on 
[a,b].  Thus  an  object  vector,  t,  can  be  constructed  as 


A 


’"’s A’r*' 


f.  =f  (i) 


(3-34) 


for  i=l,2,...,N.  Assuming  that  the  number  of  observed 
samples  is  M,  an  M-dimensional  vector,  £,  can  be  defined  as 


1 =g(i) 

i 


(3-35) 


for  i-l,2,...,M.  The  relationship  between  M and  N is  then 
given  by 


M=N+L-1 


(3-36) 


The  vector  space  formulation  of  eq.  (3-33) 


IS 


2 = Df 


(3-37) 


Where  D is  the  overdetermined  blur  matrix,  defined 


by 
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where 


(3-39) 

is  the  impulse  response  vector. 

Equation  (3-37)  states  the  fundamental  relationship 
which  exists  between  an  object  vector  and  the  cor respondinq 
observation.  Under  the  present  model  there  are  more 
observed  samples  than  unknown  parameters,  which  is  a direct 
consequence  of  the  condition  imposed  on  the  object  scene. 
This  condition  (a  scene  with  dark  setting)  is  equivalent  to 
knowledge  that  the  object  is  in  the  window.  Figure  (3-3) 
illustrates  that  equal  rate  sampling  of  the  observation  and 
the  object  results  in  more  observed  quantities  than  the 
unknown  parameters.  The  stability  of  the  system  of 
equations  defined  by  eq . (3-37)  can  be  examined  by  studying 
the  condition  number  of  D.  This  number  depends  on  the 
shape  and  the  variance  of  the  blur  function  as  well  as  the 
choice  of  the  q.f.  coefficients  c.  . Reference  [3-8] 
contains  a study  of  the  condition  number  of  overdetermined 
systems  versus  the  degrading  function  h. 

3.  5 The  Underdetermined  Model 

Continuing  on  the  discussion  of  the  previous  section, 
a more  realistic  model  evolves  if  no  restrictions  are 
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imposed  upon  the  background  of  the  scene.  The  model 
developed  in  section  (3.4)  is  rarely  very  accurate  because 
most  scenes  are  not  necessarily  placed  in  a setting  which 
has  zero  intensity,  or  even  satisfy  the  less  restrictive 
condition  of  constant  intensity.  Often,  in  the  process  of 
restoring  an  object,  the  observed  image  is  partitioned  into 
many  smaller  portions,  and  these  small  sections  are 
processed  separately.  This  partitioning  process,  by 
itself,  contradicts  any  assumption  made  on  the  setting 
which  surrounds  the  image  because  different  portions  of  an 

image  do  not  share  the  same  background  as  the  whole  image 
does. 


Assuming  lack  of  information  about  the  background  of  a 
scene,  the  object  is  assumed  to  extend  very  far  in  both 
directions,  and  consequently,  so  does  the  image.  But, 
being  able  to  handle  only  finite  segments,  one  should  be 
able  to  construct  a model  which  relates  portions  of  the 
object  to  corresponding  segments  of  the  observation. 
Figure  (3-4)  illustrates  the  concept. 

Since  the  physical  extent  of  the  observation  g(x)  is 
smaller  than  the  physical  extent  of  the  object  f (x) , the 
discrete  observation  £ is  represented  by  a smaller  number 
of  samples  than  the  discrete  object  f.  This,  of  course, 
holds  true  when  the  sampling  rate  is  kept  the  same  for  both 
the  object  and  tne  observation.  Thus,  for  an  equal 

33 


f*#  <****^'  ■ *^ 


sampling  rate,  the  system  defining  the  degradation 
phenomenon  is  equivalent  to  a set  of  linear  equations  which 
contains  a greater  number  of  unknown  parameters  than  the 
number  of  equations  available  for  solving  for  these 
parameters.  The  i-th  equation  is  given  by 


g(x.  )=  > c.  f (s)h(x  -s) 
I j i 


(3-40) 


where  x.  is  in  (a+L/2,b-L/2)  . Following  the  approach 
adopted  in  section  (3-4),  an  N-dimensional  object  vector  f, 
and  an  M-dimensional  vector  £ are  defined.  Using  notation 
B for  the  blur  matrix,  the  model  is  expressed  as 


(3-41) 


where  B is  an  MxN  matrix,  and 


M*N-L+1 


(3-42) 


Since  M<  N,  the  blur  matrix  is  not  of  full  column  rank,  and 
for  this  reason  B is  called  an  underdetermined  matrix.  The 
blur  matrix  B is  defined  as 


0 


0 


B= 


0 


^L-1  * 

L * 


• 1 


c 


0 


(3-43) 


0 c jh  L* 


Cihi 


0 

Cih^ 


Unlike  the  blur  matrix  D of  the  previous  section,  B does 
not  possess  a finite  condition  number.  Thus,  as  will  be 
discussed  in  the  next  chapter,  the  degradations  introduced 
under  the  model  of  eq.  (3-41)  are  impossible  to  remove 
completely.  Another  major  difference  between  the  two 
models  is  that  although  eq.  (3-41)  is  a more  realistic 
model  for  image  degradation  phenomena,  it  does  not  possess 
a structure  which  leads  to  the  computational  simplicity  of 
the  overdetermined  model  of  eo.  (3-37)  for  purposes  of 
image  restoration. 
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4.  NOISE  FREE  RESTORATION 


For  the  sake  of  simplicity,  let  H represent  a general 
blur  matrix  of  size  M by  N.  The  mathematical  expression 
governing  the  discrete  degradation  phenomenon,  with  H as 
the  degrading  matrix,  is  given  by 


g=Hf 


M-1) 


The  above  expression  is  a vector  space  equality,  and  in 
essence,  it  is  a set  of  M linear  equations  with  N unknown 
parameters.  The  most  straightforward  approach  for 
recovering  f,  is  to  effectively  invert  H.  If  H is  square 
and  nonsingular,  the  estimate  f is  obtained  as 


" -1 
f=H  g 


(4-2) 


where  H'^is  the  inverse  of  H.  Unfortunately,  H is  seldom  a 
square  matrix,  and  even  if  so,  H may  not  be  invertible.  To 
define  an  inverse  process  for  the  system  eq.  (4-1)  which 
will  work  in  all  circumstances,  a different  inverse  for  H 
is  defined  as 


+ T 2 -1  T 

H = lim  (H  H +d  I)  H 
d-»  0 


(4-3a) 


V ' 

• 1. 

% 

A 

% 


for  rank  N or 
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f4-3b) 


T T 2 -1 

H = lim  H (HH  +d  I) 
d-^0 


:or  rank  M,  where  H 


is  called  the  pseudoinverse  of  H 

anl  T'ae^ote.  ..e  iaentitv  It  Has  .een 

+ lA  always  exists  and 

H -^s  defined  by  eq.  (4-3),  aiway 
proven  that  H,  as  aei.  . h is 

anlcoe  ,4-.,.  THe  .ost  att.actWe  p.ope.tv  « « ^ 

stated  as  follows.  Given  the  observation  a.  the  es  « 


A + 

f=H  g 


(4-4) 


n.nrt  4-hnsp  which  minimize 
S the  vector  of  minimum  norm  among  those 


(4-5) 


i=  3”!ii 


The  above  property  of  H Introduces  the  hey  “ J"' 

.estoration  techniques  adopted  for  many  i^age  restoration 

cnprific  structure, 
Tf  H Dossesses  some  specint, 
applications.  If  H POSses 

eomputation  of  may  be  simplified  greatly,  .or  exam 
it  „ is  a square  and  singular  matrix  whose  first  r diago 

entries  are  nonzero 
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row  rank  (two 


If  H is  of  full  column  rank  or  f>. 
trequentlv  encountered  sUuations  in  image  restoration 
tasks) , simple  forms  for  H can  be  obtained.  The  next 
section  deals  with  the  problem  when  the  tank  of  a matrix  is 
equal  to  its  row  size. 


4.  1 Blur  Matrix  with  Full  Row  Rank 

Let  B denote  the  blur  matrix  possessing  full  row  rank 
Then,  the  image  degradation  model  is  described  by 


g=Bf 

The  minimum  norm  estimate  of 
eq.  (3-41)  is  given  as 

A + 

i 

where  B^  can  be  computed  from 

+ T T 2 -1 

B = lim  B (^  +d  I) 
d— *0 

Since  B has  full  row  rank, 
limit  on  the  right  hand 
out  yielding 


(3-41) 

an  object  blurred  under 

(4-8) 

(4-9) 
T 

BB  is  nonsingular , and  the 
side  of  ec.  (4-9)  can  be  carried 
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+ T T 2 t 

B =B  { lim  (BB  +d  I) } 


or  [4-3] 


+ T T -1 
B =B  ) 


l4-10b) 


Although  eq.  (4-10b)  suggests  a much  simpler  method  for 
computing  B^  than  eq.  (4-9),  an  N by  N matrix  ^ 

still  be  inverted. 

There  are  two  drawbacks  to  the  estimation  method 
described  in  this  section.  The  first  stems  from  the  fact 
that  the  model  does  not  allow  full  recovery  of  the  object 
vector  f since  there  are  fewer  equations  than  unknown 
parameters  in  the  system  of  eq.  (3-41).  The  second 

drawback  is  in  the  need  for  inverting  a matrix  of  size  N by 
N.  For  moderate  sizes  of  the  object,  the  ill  conditioning 
of  the  matrix  BB  could  cause  difficulties  [4-4].  When 
relatively  large  images  are  to  be  processed,  the  limited 
size  of  available  computers  could  put  an  intolerable 

restriction  on  the  size  of  the  object.  Often  in  situations 
like  this,  the  observation  must  be  broken  into  smaller 
segments,  each  of  which  is  used  to  estimate  the 

corresponding  object  section.  Notice  that  although  the 

A 

estimate,  f,  is  not  in  general  equal  to  the  object,  f,  the 
following  equality  always  holds 
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(4-11) 


4.  2 Blue  Matrix  with  Full  column  Rank 

Equation  (3-37)  describes  the  model,  where  D 
overdetermined  blur  matrix  of  form  given  by  eq.  (3-38) 


the 


g=Df 


(3-37) 


Since  D has  full  column  rank,  d'^D 
using  eq.  (4-3a),  the  pseudoinverse 


is  always  invertible, 
of  D can  be  obtained  as 


[4-31 


+ T -1  T 
D *(D  D)  D 


(4-12) 


Thus , 


the  minimum  norm  estimate 


is  given  by 


+ 

f = D £ 


(4-13) 


Since  system  eg.  (4-9)  has  more 
unknowns,  the  object  f can 


observed  parameters  than 
be  recovered  with  no  error. 


Thus , 


a 


(4-14a) 


or 
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(4-14b) 


A T -1  T 
f=(D  D)  D Df 


or 


f=f  (4-15) 

The  above  equality  states  the  basic  characteristic  of 
overdetermined  systems.  Full  recovery  of  the  object  is  an 
advantage  which  only  systems  of  full  column  raPK  enjoy. 
Another  superiority  of  these  systems  is  in  the  possiblity 
of  introducing  efficient  methods  which  drastically  reduce 
the  computational  complexity  of  the  associated  filters. 
The  next  section  introduces  a technique  in  which  the 
overdetermined  model  is  modified  to  pave  the  road  for 
constructing  computationally  simple  filters. 

4.  3 A Circulant  Model 

The  objective  in  this  section  is  to  establish  an  image 
degradation  vector  space  model  equivalent  to  the  one  stated 
by  eq.(3-37),  in  which  the  blur  matrix  D is  replaced  by  a 
circulant  blur  matrix  C.  A circulant  matrix  can  best  be 
explained  by  illustration:  a K by  K circulant  matrix  has 
the  following  particular  structure. 
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c c 
1 2 


• ‘k-: 


(4-16  ) 


Each  row,  or  column,  of  the  above  matrix  is  a circular 
right  shift  of  the  row,  or  column,  immediately  preceedinq. 
This  property  extends  from  the  last  row  (column)  to  the 
first  row  (column)  since  the  first  row  is  a circulant  shift 
of  the  j.ast  row  (column)  . 

To  set  up  a degradation  model  with  a blur  matrix  of 
structure  defined  by  eq.  (4-15),  two  auxiliary  vectors  are 
defined  as  follows.  Let  K be  an  integer,  where  K>M  and 
define  an  extended  object  vector  of  size  K,  f , where 

“C 


f (i)=f(i) 


for  i=l,2,..  ,N 


for  i=N+l,..  ,K 


(4-17) 


Likewise  form  an  extended  observation  vector  g of  size  K 


46 


(4-18) 


g {i)=2(i) 


for  1— »M 


=0 


for  i'-M+l , . . , K 


Next,  placing  the  impulae  response  vector  in  the  first 
column,  construct  a circulant  matrix 


h(l)  0 


h{L-l)  . h{2) 


0 h{L-l) 


C = 


h(L) 


h(L)  . 


h{l) 


(4-19) 


The  discrete  convolution  summation  of  eq.  (3-30)  is 
equivalent  to  the  vector  space  equality  [4-5] 


£ = Cf 
c c 


(4-20) 


Equation  (4-20)  which  presents  the  desired  model,  is 

(3-37),  where  D,  and  f are  replaced  by 
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similar  to  eq. 


q , C , and  f . 

“*  c 

It  is  known  that  a K by  K matrix  whi^h  possesses  K 
independent  eigenvectors  can  be  diagonalized  through  a 
similarity  transformation  [4-6]  . Hunt  [4-7]  demonstrates 
that  matrices  of  the  form  of  eq.  (4-15)  have  K independent 
eigenvectors,  and  that  the  Fourier  transform  diagonalizes 
circulant  matrices.  In  fact,  the  eigenvectors  of 
circulants  are  the  Fourier  basis  vectors.  Let  A represent 
the  Fourier  transform  matrix,  thus 


2rij 

A(i,k)=exp{-(— ik] 


(4-21) 


The  similarity  transform  which  diagonalizes  the  blur  matrix 
C of  eq.  (4-19)  is 


C=S'^A  A 


(4-22) 


where  A is  the  diagonal  matrix  of  the  eigenvalues  of  C 


(4-23) 


The  X are  obtained  by  Fourier  transforming  the  first 
i 

column  of  C.  Thus, 
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(4-24) 


L-  1 


2nj 


=y  ^exp{-(i-l)k }h(k+l) 

K 


for  i— 1,,,.,K.  To  study  SQ.  (4—20)  in  Fourier 
substitute  eq.  (4-22)  in  eq.  (4-20)  to  obtain 


g =A  A Af 
*c  — - — c 


Next,  eq.  (4-27)  is  rearranged  as  the  following 


K.. 


M = AAf 
— ^ 


By  definition  ^ and  ^ are  the  discrete 

transforms  of  the  vectors  g and  f , respectively. 

c — c 


F =Af 
~ C c 


G =Ag 

~ C — 


where 


2nJ 


= _(k+l)exp{-(i-l)  k-^} 

k=0 
K-1 


for  i=l,...,K.  Thus,  in  the  transform  domain. 


eq. 


space , 
(4-25) 

(4-26) 

Fourier 

(4-27a) 

(4-27b) 

(4-28) 

(4-20) 
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r 


simplifies  to 


(4-29) 


G^= 


Since  A is  a diagonal  matrix,  eq.  (4-29,  is  actually  a 
scalar  equation 


G^(i)  = X Fe^i) 


(4-30) 


p multiply  both  sides  of 
for  i=l,...rK.  TO  restore  , multiply 

eq.  (4-29)  by  A^.  Thus, 


A -1 
F » AG 
-c  — c 


(4-31) 


or  in  the  scalar  form 


F (i)=--G  (i) 
c A . c 

^ i 


(4-32) 


Inverse  Fourier 


A ,A 


i*l r . . »K 

transforming  of  F^  results  in  the  estimate 

(4-33) 


...  . f can  be  obtained  by  extracting  out 

The  object  estimate,  £,  can  e 


the 


first  N entries  of  fc* 


i=iiS  ie 


(4-34) 
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where  is  a selection  matrix  of  the  following  form 


K 


I I 0 


N 


I 


(4-35) 


Figure  (4-1)  illustrates  the  relationship  between  the 
models  expressed  by  eq.  (3-38)  and  eq.  (4-20). 

As  noted  earlier,  the  model  developed  in  this  section 

is  equivalent  to  the  overdetermined  model  established  in 

section  (3.  4);  thus,  the  restriction  of  objects  with  dark 

background  still  exists.  The  incentive  for  taking  a detour 

by  introducing  the  circulant  model  is  in  the  present 

system's  computational  superiority  over  the  previous  model. 

The  inverse  filter,  a' used  in  estimating  f does  not 

~c 

require  a matrix  inversion.  Thus,  the  ill-conditioning,  or 
the  large  size  of  the  system  is  not  a major  obstacle  in  the 
process  of  computing  the  circulant  filter. 

4.  4 Experimental  Results 

Figure  (4-2)  illustrates  an  object  with  zero 


51 


background.  The  object  in  this  case  is  a set  of  bars  of 
constant  intensity  separated  from  each  other  by  a set  of 
bars  with  zero  intensity.  The  object  and  its  background 
form  a square  picture  which  is  sampled  at  256x256  uniformly 
spaced  points.  Figure  (4-3a)  shows  this  object  after 
undergoing  a one-dimensional  blur  degradation  of  Gaussian 
shape  with  standard  deviation  of  2.5  pixels.  Figure  (4-3b) 
shows  the  estimate.  This  estimate  has  been  obtained  using 
the  method  explained  in  section  4.  3.  Since  the 
restoration  has  been  performed  in  the  Fourier  domain  (a 
scalar  operation) , the  relatively  large  size  of  the  image 
is  of  no  concern.  The  estimate  is  identical  to  the  object 
itself.  This  interesting  achievement  holds  for  all 
examples  in  which  the  object  possesses  a black  background, 
and  the  environment  is  noise  free.  Figure  (4-4a) 
illustrates  the  object  after  undergoing  an  extreme  amount 
of  blur.  The  degradation  models  a motion  blur  for  which 
L*15.  Figure  (4-4b)  is  the  estimate  which  is  error  free. 


If  the  images  discussed  above  were  to  be  restored 
using  an  under  determined  model  (section  4.  1),  the  large 
size  of  the  image  would  require  that  the  observation  be 
partitioned  into  smaller  segments.  Note  that  an  estimate 
of  th-!  form  given  by  eq.  (4-8)  cannot  be  obtained  by  direct 
Fourier  techniques  [4-5],  [4-8].  Figure  (4-5)  plots  the 
error  for  the  underdetermined  model  estimate.  Each 
observed  segment  has  25  pixels  and  is  assumed  to  estimate 
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17  pixels 
different 
deviation 
to  be  9. 


of  the  object.  The  error  is  plotted  for  two 
blur  models:  a Gaussian  shaped  blur  of  standard 
2.5;  and  motion  blur.  In  both  cases  L is  assumed 
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5.  WINDOWING  OPERATOR 


The  previous  section  introduced  a circular  image 
degradation  model  which  resulted  in  a computationally 
efficient  restoration  technique.  Unfortunately,  this 
technique  can  only  be  applied  to  a scene  which  possesses  a 
dark  background.  Hence,  portions  of  a given  image,  or  an 
object  placed  in  a nonzero  intensity  setting  cannot  be 
recovered  by  brute  forcing  the  circulant  filter  on  the 
observed  image.  in  fact,  since  this  kind  of  imaging  can 
not  be  expressed  by  an  overdetermined  model,  employment  of 
a system  of  the  form  of  eq.  (4-20)  for  representing  the 
degradation  phenomenon,  when  objects  of  unknown  background 
are  involved,  can  cause  a modelling  error  at  the  boundaries 
of  the  image.  This  error  can  be  looked  at  as  signal 
correlated  noise;  hence,  considering  the  ill-conditioning 
inherent  in  imaging  models,  the  implementation  of  the 
circulant  filter  becomes  catastrophic.  To  illustrate  this 
claim,  figure  (5-la)  is  selected  as  a test  image.  This 
image  is  obtained  by  slightly  blurring  the  original  object. 
Figure  (5-lb)  shows  the  same  image  after  an  attempt  is  made 
to  restore  the  center  portion  of  the  object  by  utilizing 
the  fast  Fourier  technique  of  eq.  (4-31).  Observing  the 
high  frequency  noise  component  in  figure  (5-lb),  it  is 
evident  that  an  unwise  choice  of  a degradation  model  not 
only  does  not  improve  the  observation,  but,  on  the 
contrary,  it  could  be  quite  destructive.  To  approach  this 
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problem  in  a more  systematic  manner,  the  relationship 
between  the  background  of  the  object,  the  system  model,  and 
the  observed  image  must  be  studied.  The  next  section 
discusses  this  subject  and  suggests  means  of  treating  a 
given  observation  in  order  to  modify  the  physical  samples 
to  suit  a desired  model. 

5.  1 Overdetermined  System  with  Unrestricted  Observation 

The  primary  difference  between  the  two  sets  of 
equations  (3-41)  and  (3-37)  is  stated  as  follows.  There 
are  only  M equations  (or  observed  parameters)  in  the  system 
of  eq.  (3-41)  to  solve  for  M+L-1  unknown  parameters,  thus, 
since  L>1,  an  exact  restoration  of  the  object  is 
impossible.  On  the  other  hand,  the  system  described  by 
eq.  (3-37)  has  M observations  for  only  m-L+1  unknowns. 
This  means  that  number  of  equations  is  actually  larger  than 
the  minimum  number  necessary  for  a complete  restoration  of 
the  object.  Now,  considering  the  assumption  made  on  the 
data  of  the  latter  model,  it  appears  that  the 
overdetermined  model  is  in  essence  equivalent  to  its 
underdetermined  counterpart  if  some  of  the  unknown 
parameters  of  the  underdetermined  system  are  set  to  be 

equal  to  zero.  An  M by  n underdetermined  system  is 
represented  by 
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jjesB  w 


q=Bf 

~ (3-41) 

The  following  lemma  states  the  oompatiblity  of  system 

eg.  (3-41)  with  an  overdetermined  system  of  form 

eg.  (3-37),  where  the  corresponding  blur  matrix  is  of  size 
M by  M-L+1. 


Lemma  5-1.  The  underdetermined  system  of  eg.  (3-41) 
becomes  egnivalent  to  an  overdetermined  system  of  the  form 
eg.  (3-37)  if  the  first  and  the  last  L-1  entries  of  the 
object  vector  f are  set  egual  to  zero. 

Proof.  Let  e represent  the  object  vector  after  the  first 
and  the  last  L-1  entries  are  set  to  zero. 


e(i)*f  (i) 


if  L<i<N-L+l 


e(i)«0 


otherwise 


(5-1) 


and  assume  d,  a vector  of  size  M-L+1, 
nonzero  center  portion  of  e.  The  observati 
to  object  e is  given  as  follows 


represents  the 
on  corresponding 


3=^ 


Notice  that 


(5-2) 
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Be=Dd 


(5-3) 


where  D is  the  overdetermined  matrix  of  dimensions  M by 
M-L+1  and  is  given  by  eq.  (3-34).  Equality  (5-3)  holds 
because  of  the  particular  structure  that  both  D and  B have. 
If  eq.  (5-3)  is  substituted  in  eq.  (5-2)  then 


q=Dd 


(5-4) 


Equation  (5-4)  is  the  desired  result.  To  find  g,  subtract 
eq.  (3-41)  from  eq.  (5-2)  giving 


2-g=B(e-f ) 


(5-5) 


or 


3=2“B(f-e) 


(5-6) 


Introducing  an  N by  N selection  matrix 


(5-7) 


eq.  (5-6)  simplifies  to 


a=q-BS^f  (5-8) 

- — N- 

Equation  (5-8)  holds  true  since 

f-e=S^f  (5-9) 

N- 

The  role  of  is  to  select  the  first  and  the  last  L-1 
-N 

entries  of  f.  Vector  f-e,  in  fact,  represents  the 

background  of  the  object,  and  vector  B(^-e)  in  eq.(5-6) 
represents  the  contribution  of  this  background  to  the 
image . 

What  lemma  (5-1)  suggests  is  restated  as  follows.  Any 
observed  image  can  be  suitably  processed  for  an 

overdetermined  modei  provided  that  the  intensity  function 
describing  the  setting  of  the  object  is  obtainable.  The 
major  drawback  of  the  above  statement  is  that  the  intensity 
function  of  the  surrounding  of  the  object  is  not  usually 
known  a priori.  But  often  this  function  can  be  estimated 
with  an  acceptable  accuracy.  This  is  the  subject  of  tne 
following  section 

5.  2 Windowing  of  the  observation 

Since  the  only  available  source  of  information 
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concerning  the  scene  is  the  observed  vector , the  estimator 
which  estimates  the  object  background  must  take  use  of  the 
the  physical  sample  vector  Assume  that  matrix  W 
represents  the  combination  of  a background  estimator  plus 
the  system  which  removes  the  effect  of  the  estimated 
background  from  the  observation.  The  product  of  the  image 
vector,  by  matrix  W results  in  the  desired  observation, 
3,  which  can,  successfully,  be  used  in  an  overdetermined 
system.  The  structure  of  matrix  W (the  windowing  matrix) 
is  explained  in  the  following  paragraph. 


Let  3 represent  the  observation  vector  if  the  first 
and  the  last  L-1  entries  of  f are  zero.  The  objective  here 
is  to  express  3 in  terms  of  the  elements  of  the  observation 
3 according  to  the  relation 


q(i)  = 


g(i)-/  ^h(L+l-j)  f (j  + i-l) 

L-l-M+i 

[ g(i)  i:  h { j ) f (N+1- j + i-M) 
j = l 


if  i<L 

L<i<M-L  (5-10) 

if  i>M-L+l 


Since  the  entries  of  £ are  not  known,  the  correspondence  of 
eq.  (5-11)  cannot  be  made  directly.  However  by  making  an 
assumption  on  the  continuity  of  the  original  image  vector 


that 


f(i)=f(L)  for  i<L 


(5-11  ) 

6 5 


and 


,<**  ^ 5»», 

--^  A" 


f (N-i+l)=f (N-L)  i>N-L 


then  the  vector  c[  can  be  estimated  as 


a=w2 


(5-13) 


where  W is  an  M by  M matrix  of  the  form  given  by 
eq.  (5-14). 


A zero  order  predictor  is  inherent  in  the  structure  of 
matrix  W expressed  by  eq.  (5-14).  The  prediction  of  the 
image  background  is  the  main  idea  in  expression  (5-13). 
Therefore,  the  success  of  the  operation  defined  by  matrix  W 
depends  on  the  validity  of  the  prediction  method  used  to 
obtain  W.  There  are,  of  course,  other  prediction  methods 
which  can  be  employed.  For  example  a first  order  predictor 
results  in  a smaller  (overall)  mean  square  error.  Figure 


. I|9^.  **S  ."t 


a 


1 


0 


(5-14) 


(5-2)  illustrates  the  expected  mean  square  restoration 
error  of  the  object-estimate  for  two  prediction  alqorithms, 
as  a function  of  element  correlation  p.  Fortunately  the 
zero  order  predictor,  in  practice,  produces  sufficiently 
accurate  results.  Since  this  predictor  has  a simple  form, 
it  has  been  employed  in  the  remaining  material  of  this 
text . 
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It  is  interesting  to  note  that  the  physical 
contribution  of  multiplying  an  image  vector  by  W is  an 
attenuation  of  the  first  and  the  last  L-1  entries  of  the 


corresponding  vector.  This  is  expected  since  the 
observation  resulting  from  the  class  of  objects  encircled 
by  a dark  region  illustrates  dim  object  boundaries 
(boundary  of  the  object  with  its  background) . Figure 
(5-3a)  represents  the  general  form  of  a typical  observed 
image  line  while  figure  (5-3b)  shows  the  same  image 
function  after  undergoing  a windowing  operation. 


5.  3 Error  Analysis 


Assume  X,  a vector  of  size  N,  represents  the  object, 
and  let 


2=^ 


(5-15) 


symbolize  the  observed  image  after  x has  undergone  a 
degradation  of  known  impulse  response.  The  size  of  £ is 
given  by 


M=N-L+1 


(3-38) 


The  objective  hare  is  to  estimate  the  center  elements  of  x 


using  the  overdetermined  system  model.  Since  the  physical 
sample  vector  £ and  the  system  equation  (3-41)  are  not 
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co.7.patible,  the  observed  vector  3 must  undergo  the 
windowing  operation  described  in  the  previous  section.  The 
modified  observation,  q,  is  obtained  as 


a=W2 


(5-16) 


or 


g=WBx 


(5-17) 


At  this  stage,  the  filter  derived  for  the  overdetermined 

system.  See.  4-2,  can  be  employed  to  estimate  the  center 
part  of  X.  Tiierefore, 


A + 
f=D  q 


(5-18) 


Where , 
center 


D*is  given  by  eq.  (4-12)  , and  f Is  the  estln,ated 
part  of  X.  The  length  of  f , k.  Is  given  by 


K=M-L+1 


(5- .19) 


If  expression  (5-17)  is  used  to  substitute  for  g 
eq.  (5-18)  becomes 


(5-2f) 


A -f 
f=D  WBx 


f,  a vector  of 
image  vector 
"Sing  the  selection 


size  K,  denote  th 

— * This  Vector 
■^3tr  ix 


® center  portion  of  th 
be  extracted  from 


e 

X 


K 


I 


(5-21) 


Thus, 


;:=S2^  X 
- -=N  ^ 


(5-22) 


figure  (5-41  ,•  1 1 

4)  Illustrates  the  oo 

'’cctor,;  X,  2,  3,  and  (.  ■^'^^POndence  between  the 


The  eeti„at;on  error,  e ie  a.,- 

' - <Jei'ned  as 


A 

e=f-f 

(5-23) 

hhe  error  vector  fro.  , ,,,,, 

the  error  vector  e la  =‘“tetrcal  point  of 

assumed  to  ho  m,. 

be  mean  zero,  this 
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Figure  (5-4)  The  relationship  between  the  physical  observed  samples 
and  the  object. 


assumption  is  merely  made  for  the  sake  of  convenience.  It 

is  always  possible  to  subtract  the  mean  value  of  the  object 

from  X.  This  will  automatically  modify  all  the 

corresponding  vectors  to  become  mean  zero.  The  error 

covariance  matrix,  R , is  defined  by 

-e 

R =E{ee'^}  (5-24) 

~e  — 

Where,  the  notation  E denotes  the  ensemble  average.  Using 

the  expression  for  e from  eq.  (5-23),  R becomes 

“ — e 

R =E{f-f  ) (f-f  )^}  (5-25) 

— e ~ ^ ~ 


or 


(5-26) 


Let 


R =E{xx^}  (5-27) 

~x  — 

represent  the  correlation  matrix  of  the  object.  The  right 
hand  side  of  eq.  (5-26)  can  then  be  evaluated  element  by 
element  as 

E{ff'’^}=^^  E{^'’^}  (^^  )'^  (5-28) 
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E{ff  }=S2  R (S2  ) 

~ — N X — X 


The  second  term  can  be  evaluated  as 


aaT  + rj.  T T +T 

E{^  }}=D  TO[E{)«  }]B  W (D  ) 


or 


aaT  + T T ^ T 
E{j[f  }^D  WBR  B W (D  ) 


T 

The  expected  value  of  ff  can  be  obtained  as 


AT  + 'r  K 

E{^  }»'D  ^[E{^  H{S2^.  ) 


or 


at  + r 

E{ff  }=D  WBR  ) 

^ N 


aT 

Likewise  E{^  } can  be  derived  as 


A T K T T + T 
E{ff  }=S2*,  R B W (D  ) 
— — X ' 


Now  the  total  error  covariance  can  be  expressed  as 


TTTK  KT+^ 

R^=D^WBRyB  W )”D  WBR 


N 


(S2^ 
— N 


-S2j^B^W^(D+) 


(5-35) 


It  is  possible  to  process  the  physical  samples  of  the 
blurred  image,  g,  with  the  filter  derived  under  the 
underdetermined  system  model  assumption.  In  this  case  the 
alter  is  given  by  eg.  (3-41).  Therefore,  the  estimate,  x, 

is  given  by 


2 


(5-36) 


A 

And  the  center  portion  of  the  estimate,  f,  is  simply 
obtained  by  premultiplying  eq.  (5-36)  by  to  yield 


A ka 


(5-37) 


or 


(5-38) 


The  error  term 


A 

e=f-f 


(5-39) 


has  the  following  covariance 
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A at 
R =E{ (f-f) (f-f)} 


(5-40) 


Going  through  the  steps  similar  to  those  described  in  the 


previous  case,  the  error  covariance  can  be  derived  as 


Iv  ^ 


(5-41) 


where,  I is  the  M by  M identity  matrix.  Figure  (5-2)  of 
the  previous  section  illustrates  the  expected  mean  square 


restoration  error  of  f as  a function  of  the  correlation  of 


elements  in  f under  the  assumption  that  f is  a sample  of  a 
Markov  process  with  correlation  factor  p. 


Figu..e  (5-2)  illustrates  that,  except  when  the  object 


element  correlation  coefficient  is  near  unity,  the  error 


resulting  from  the  estimate  given  by  eq.  (5-20)  is  larger 


than  the  one  resulting  from  eq.  (5-38).  But,  consideLing 


that  most  images  observed  by  a human  viewer  possess  strong 


correlation  among  their  sampled  pixels,  this  extra 


contribution  of  error  is  not,  usually,  unreasonably  high  in 


practice.  Also,  with  the  D operator,  contained  in 
eq.  (5-20),  it  is  possible  to  perform  the  restoration  by 
Fourier  domain  processing  quite  efficiently. 
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5.  4 Experimental  results 


Figure  (5-5a)  illustrates  an  image  corrupted  by  a 
Gaussian  blur  modelling  a blur  degradation  caused  by 
imaging  through  a turbulent  atmosphere.  The  center  portion 
of  the  image  has  been  filtered  to  produce  the  corresponding 
section  of  the  object.  A zero  order  predictor  was  employed 
for  estimating  the  object  background,  and  the  restoration 
is  performed  in  the  Fourier  domain  as  indicated  in 
Sec.  5.  3.  Figure  (5-5b)  shows  the  image  after  the  central 
portion  of  the  object  hcs  been  restored.  Unfortunately, 
since  the  background  predictor  is  not  error  free,  the 
i^ilter  is  not  applicable  for  severe  amounts  of  blur.  Since 
the  system  modelling  the  degradation  process  is  basically 
ill-conditioned,  the  restoration  technique  greatly 
amplifies  any  uncertainty  in  the  observation.  Figure 
(5-6a)  illustrates  a test  object  after  undergoing  a severe 
amount  of  blur.  Figure  (5-6b)  is  an  attempt  to  restore  the 
object,  which  clearly  has  been  unsuccessful.  If  the 
background  of  the  object  is  of  constant  intensity,  the 
zero-order  predictor  inherent  in  the  windowing  matrix  can 
make  an  exact  estimate  of  the  background.  This  would 
insure  an  error  free  observation.  Figure  (5-7a) 
illustrates  this  case.  The  center  part  of  the  object  has 
been  processed  to  artificially  generate  a constant 
intensity  background.  Figure  (5-7b)  is  the  restored 
version  of  fig.  (5-7a).  As  before,  only  the  center  part  of 
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6.  NOISY  RESTORATION 


The  models  formulating  the  image  degradation 
phenomenon  in  the  preceeding  chapters  have  ignored  the 
presence  of  noise.  In  practice,  however,  imaging  systems 
are  seldom  totally  noise  free.  In  image  forming  systems, 
noise  or  uncertainty  may  arise  from  a variety  of  sources; 
probably  the  most  common  sources  of  noise  are  measurement 
and  recording  errors.  Scanners,  the  devices  which  measure 
images,  invariably  add  an  element  of  uncertainty  to  the 
measured  (or  scanned)  signal  [6-1].  Usually,  after  an 
image  is  scanned,  it  is  operated  upon  in  a computer  of 
finite  precision,  creating  truncation  errors.  Coding  and 
channel  errors  are  caused  if  the  particular  image  is  to  be 
transmitted  through  a noisy  channel  [6-1).  Lastly,  film 
noise  m<jy  be  added  to  the  image  when  the  signal  is  recorded 
[6-2].  It  is,  of  course,  impractical  to  make  an  exhaustiv^e 
list  of  all  noise  producing  elements  in  an  image  forming 
process,  but  the  noise  sources  listed  above  are  the  most 
significant. 

The  next  section  studies  the  continuous  image 
degradation  problem  in  the  presence  of  additive  white 
noise . 
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6.  1 The  Continuous  Model 


The  continuous  image  deqra<  ^.tion  model  of  eq.  (2-1)  is 
modified  in  the  presence  of  additive  white  noise  to 


f(s)h(y-s)ds  +n(x) 


(6-1) 


where  n(x)  represents  the  noise  component.  As  a result  of 
the  existence  of  noise  in  the  system,  the  inverse  filtering 
technique  cannot  be  employed  to  recover  the  object.  If  any 
attempt  is  made  to  force  the  inverse  filter  to  process  the 
physical  sample  image  g(x),  the  high  frequency  component  of 
the  noise  will  be  greatly  amplified  ruinin’  the 
restoration.  To  avoid  high  frequency  amplification  of  the 
noise,  the  inverse  filter  can  be  truncated  at  a certain 
point.  The  point  can  be  selected  so  that  beyond  this  point 
the  noise  power  exceeds  the  signal  energy  [6-3]  . Figure 
(6-1)  illustrates  tnis  metliod. 

A more  intelligent  technique  to  recover  the  object  in 
the  presence  of  noise  is  the  classical  Wiener  filter.  This 
filter  controls  the  noise  component  by  keeping  into  account 
the  ratio  between  the  signal  power  and  the  noise  variance 
at  each  point  of  the  Fourier  space.  Tfie  Wiener  filter 
output  is  a minimum  mean-square  error  estimate  of  the 
original  object,  provided  that  the  statistics  used  in  the 
filter  ai.^  carefully  obtained. 
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6.  2 The  Discrete  Wiener  Filter 


In  the  presence  of  noise,  the  system  governing  the 
image  degradation  phenomenon  is  formulated  as  16-31 


£=Hftn 


(6-2) 


where  H is  a blur  matrix  and  n is  the  noise  component 
vector.  TO  derive  the  minimum  mean-square  error  filter  for 
the  system  described  by  eq.  (6-2)  the  «ell  known 
orthogonality  principle  (6-4| , (5-61  can  be  employed. 
Letting  U represent  the  desired  filter;  the  estimate  f is 

obtained  as 


A 


(6-3) 


According  to  the  orthogonality  principle  the  following 
equality  must  hold 


^ T „ 
E(f-f)£  =0 


(6-4) 


After  substituting  for  f from  eg.  (6-3)  and  carrying  out 
the  appropriate  steps,  the  solution  is  given  as  [6-3] 


U=R  H^(H^  H^+V)"^ 
f - — f - 


(6-5) 


;re  R is  the  correlation  matrix  of  the  object  and  V is 


the  noise  correlation  matrix.  The  object  and  the  noise  are 
assumed  to  be  uncorrelated.  The  following  matrix  identity 
[6-6] 


( c+^e'^r  ^ ( A c ’b  j ^ B^C"  ^ 


(6-6) 


can  be  used  to  obtain  a different  formulation  for  filter  U 
given  by 


-1  -1  -1  T -1 

U=(R  +H^ V H)  H T 

- — f _ _ _ - - 


(6-7) 


The  error  variance  matrix  C is  given  as  [6-7] 

— e 


T T T -1 

t ^ \ ftl 


C =^-(R  hM  (HR  H^+V)  (R,H  ) 


(6-8) 


And  by  using  the  matrix  identity  [6-6] 


-1T“1  T T”1 

(C  +B  ^)  =C-CB  (B£B  +A)  BC 


(6-9) 


the  error  term  simplifies  to 


£e  = ( R ^+H^V  ^)'  ^ 


(6-10) 


For  a white  noise  process,  the  noise  correlation 
matrix  has  the  following  simple  form 
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V=d^  I 


(6-11) 


where  d ^ is  the  noise  power,  and  I denotes  the  identity 
matrix.  Representing  the  image  power  by  i , expression 
(6-7)  becomes 


h'^H) 


1 T^-2 
H d 


or 


-1  T -1  T 
U=(S  C +H  H)  H 


(6-13) 


where  represents  the  normalized  object  covariance 

matrix,  and  S represents  the  signal  to  noise  ratio 

_ ^ (6-14) 

Likewi-e,  expression  (6-5)  simplifies  to 

U=C  h'^(HQ,  h’^+S  ^I)'^  (6-15; 

f - — f - - 


Although  expressions  (6-13)  and  (6-15)  are  equivalent,  the 
first  expression  is  employed  when  blur  matrix  H is  of  full 
column  rank  (type  D) , and  the  second  one  is  used  if  H has 
full  row  rank  (type  B)  . The  reason  for  this  is  simply 
because  of  computational  savings.  Note  that  for  an 
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overdetermined  matrix  the  number  of  columns  is  less  than 

Lhe  number  of  rows  (Sec.  3.4).  Hence,  the  dimensionality 

T 

of  (NxN)  is  less  than  the  dimensionality  of  ^ (MxM)  . 

And,  for  an  underdeterir  ’ ned  matrix,  the  order  is  reversed; 
the  dimensionality  of.  (NxN)  is  less  than  the 

dimensionality  of  (MxM).  Thus,  for  the  case  of 

overdetermined  matrices  an  N by  N matrix  (N<M)  has  to  be 

inverted,  and  in  the  other  case  (underdetermined  blur 
matrix)  the  inversion  of  an  M by  M matrix  (M<N)  is 
required . 

The  necessity  of  inverting  matrices  of  relatively 
large  size  is  an  unattractive  property  of  the  discrete 
Wiener  filter.  To  prevent  this,  the  next  section 

introduces  a technique  which  eliminates  the  matrix 
inversion  requirement. 


6.  3 The  Fast  Wiener  Filter 

Section  (4.  3)  introduced  an  image  degradation  model 
which,  by  utilizing  the  Fourier  domain  properties  of 
circulant  matrices,  successfully  eliminated  the  matrix 
inversion  requirement  of  the  pseudoinverse  filter. 
Although  the  circular  model  provides  attractive 
computational  savings  in  a noise-free  environment,  iu 
cannot  be  applied  to  a noisy  observation.  In  fact,  since 
the  windowing  operation  of  Sec.  (5.2)  introduces  an  element 
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I 

t 
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of  uncertainty  in  the  observation  space,  the 
pseudoinverse  technique  i&  unable  to  recover  the  object 
after  it  has  undergone  a severe  amount  of  degradation.  To 
overcome  this  shortcoming,  it  is  necessary  to  consider  the 
existence  of  noise  in  the  system.  Thus,  this  section  is 
devoteu  to  the  derivation  of  a system  model  which  retains 
computational  simplicity  and,  at  the  same  time,  tolerates 
noise-corrupted  observation  samples. 

In  the  presence  of  additive  noise,  the  image 
degradation  model  with  a circulant  blur  matrix  C is 
formulated  as 


f. 


4 


2 =Cf  +n  (6-16) 

c c 

According  to  the  developments  and  the  definitions 

established  in  Sec.  (4.  3),  both  vectors  £ and  f are  of 

size  K (K>M) , and  represent  the  hypothetical  observation 

and  object,  respectively.  Since  the  actual  physical  sample 

image,  £,  is  of  size  M and  the  actual  object  vector,  f,  is 

v.f  size  N,  the  last  K-M  entries  of  g are  pure  noise,  and 

c 

the  last  K-N  entries  of  fc  are  deterministic  and  equal  to 
zero.  In  searching  for  a minimum  mean-square  error  filter 
Corresponding  to  the  observation  in  eq.  (6-16),  the  Wiener 
filter  may  seem  the  logical  candidate;  however,  this  is  not 
true.  To  illustrate  this  fact,  consider  the  Wiener  filter 
solution  for  the  white  noise  case  as  given  by 
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-1-1  T - 1 T 
U=(S  C +C  C)  C 


(6-17) 


where  C is  the  covariance  matrix  of  the  object.  Since  the 


covariance  matrix  C,  is  not  circulant,  the  undesirable  KxK 


matrix  inverse  operation  cannot  be  avoided  by  simply 


Fourier  transforming  the  filter.  Furthermore,  since  some 


of  the  entries  of  ^ are  deterministic,  C is  a singular 

1 ^ 

matrix,  and  thus  does  not  exist.  The  preceeding 


argument  implies  that  a circulant  Wiener  filter  is  not 


achieveable  because  of  the  particular  statistical 


properties  of  the  vector  f . It  is  conceivable  that  a 


circulant  and  nonsingular  covariance  matrix  could  be 


obtained  if  the  postulated  object  vector,  f^  , had  a proper 


statistical  background.  To  obtain  a nonsingular  covariance 


matrix,  none  of  the  entries  of  f^  can  be  deterministic.  To 


achieve  this  result,  form  a new  vector  ^ by  augmenting  the 


actual  object  f with  a vector  ^ which  has  certain  desired 


statistical  properties  to  be  described  later.  The 


augmented  vector  is 


(6-18) 


The  covariance  matrix  of  the  above  vector  i: 


f 


T 

E(f  f )= 
-c  — c 


C,  Cf 
-f  -fy 


C , C 
-yf  -y 


where 


C =E(ff  ) 


is  the  actual  object  NxN  covariance  matrix,  and 


=E(fx  ) 

— fy  — 


C =E(Yy^)  (6-22) 

~y 


is  the  (K-N)x(K-N)  dimensional  covariance  matrix  of  At 

this  stage,  the  aim  is  to  illustrate  that  a certain 


statistical  assumption  on  ^ plus  a specific  value  of  K can 

result  in  a circulant  and  positive  definite  matrix  C,. 

'c 


As  an  example  let  the  object  vector,  f,  consist  of 

I 

exactly  four  samples  given  by 


Assuming  that  ^ arises  from  a Markovian  random  process,  the 
covariance  matrix  of  f has  the  following  form 


1 p 

P 1 P 

2 

p p 1 p 

p p p 1 


(6-24) 


where  p is  the  element  correlation  coefficient.  Since  the 
size  of  f,  N,  is  four,  is  of  size  4x4.  Selecting  K to 
be  equal  to  6 , a 6 by  6 circulant  matrix  is  introduced 


1 P P^  P^  P^  P 

P 1 P E?  P^  P^ 

P^  P 1 P P^  P^ 

p^  F?  p 1 p p2 

P^  P^'  F?  P 1 P 

P P^  P^  P^  P 1 


(6-25) 


Since  the  matrix  in  eq.  (6-25)  is  symmetric,  it  must  be 

proven  that  it  is  positive  definite  in  order  to  conclude 

that  C is  indeed  a covariance  matrix.  Since  C is 
~fc  — f 

circulant,  Fourier  transformation  of  the  first  column  can 
generate  the  eigenvalues  of  this  matrix  [6-8].  Figure 
(6-2)  lists  all  the  eigenvalues  of  according  to  their 
number.  Since  the  element  correlation  is  always  less  than 


'4 


f 


/ 
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unity,  all  the  eigenvalues  are  positive.  Thus,  C^.  is  a 


positive  definite  and  circulant  matrix.  Since  K=6  in  this 


example,  vector  ^ is  exactly  of  size  2 


l=[Vyy2^ 


(6-26) 


According  to  eq.  (6-25),  ^ must  have  the  following 

statistics  (^  is  mean-zero) 


1 P 


P 1 


(6-27) 


p2  P 


3 2 
P P 
2 3 

P P 

2 

P P 


(6-28) 


Equations  (6-27)  and  (6-28)  can  be  used  to  generate  the 


random  process  samples 


To  generalize  the  above  example,  the  case  of  an 


N-d imensional  object  is  considered  here.  Parameter  K is 


selected  according  to  the  following  rule 


K=N+N-2 


(6-29) 


And  similar  to  the  covariance  matrix  of  eq.  (6-25) , a KxK 
covariance  matrix  is  constructed  using  the  NxN  Markovian 
matrix.  To  acheive  this,  N-2  center  entries  of  the  first 
row  are  folded  over  and  attached  to  the  first  row  itself  to 
produce  the  first  line  of  the  KxK  circulant  covariance 
matrix.  Since  the  resulting  matrix  is  circulant,  only  the 
first  row  is  needed  to  construct  the  remainder  of  the 
matrix.  Figure  (6-3a)  illustrates  the  first  row  of  an  NxN 
Markovian  matrix,  and  fig.  (6-3b)  shows  the  resulting  first 
line  of  the  circulant  matrix.  Figurt  (6-4)  illustrates  the 
resulting  KxK  circulant  matrix.  The  eigenvalues  of  this 
matrix  can  be  obtained  by  Fourier  transforming  the  first 
row  of  the  matrix  itself. 

Let  a K-dimensional  vector  £,  which  represents  the 
first  row  of  the  circulant  matrix,  be  defined  as 


r=[  1 


p2  ...  pN-2..  p2  Pi 


(6-30) 


The  Fourier  transform  of  £,  R,  is  a vector  of  size  K which 
is  obtained  from 


R=Ar 


(6-31) 


where  A represents  the  discrete  Fourier  transform  matrix. 
The  vector  R contains  the  eigenvalues  of  the  matrix  in 
fig.  (6-4).  Since  £ has  a particular  structure,  it  is 
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? ... 


First  line  of  an  N by  M Markovian  covariance  matrix 


- , <u 

( 1 P P 


. g-3  pN-2pN.l  ^-2pN-3__  p2  p j 


First  row  of  the  resulting  Markovian  circulant  matrix 


Figure  (6-3)  The  first  line  of  an  NxN  Markovian  matrix 
giveo  rise  to  the  first  line  of  a circulant  matrix. 
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possible  to  find  d closed  form  solution  tor  the  entries  of 
R.  By  definition  of  the  Fourier  transform,  the  (k+l)th 
entry  of  R,  R(k+1) , is  given  as 
2N-3 


R(k*l).2j^(jtl)expt-l-j^^nkl 


(6-32) 


j=0 


where  k varies  from  zero  to  K=2N-3,  and  i is  the  square 
root  of  unity.  Equation  (6-31)  can  be  partitioned  as 
N-1 

« > 2ni  1 • 1, 1 ■ 

R(k+1)=^  t (j  + l)exp(-[-j2Nr^l  3k)  + 


zni 


r(j+l)exp{-[-(2Nr^lik} 


(6-33) 


The  exact  value  of  r(j+l)  can  be  obtained  from  eq.  (6-30). 
If  eq.  (6-30)  is  used  in  eq.  (6-33),  then  along  with 
further  simplifications  the  following  equation  can  be 

obtained 


k N-  i 

R(k+l)=l+ (-1)  P 


N-2 

+ ^ V {expl- 

j = l 


2N-2  2N-2 


(6-34) 


or 


N-2 

k N-1  j _ , 2njk  ^ 

R(k+l)  = l+(-l)  p +2^  V Cos[ — ] 


2N-2 


(6-35) 
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Obcer'^ing  the  following  identity  [6—9] 


N+l 


ll-pCos(t)]  [1-p  Cos(Nt)  ]+p  Sin(t)Sin(Nt) 


l-2pCos 


(t)+f? 


(6-36) 


l-pCos(t)+p^+Cos(N-l) t-p  Cos(Nt) 

2 

l-2pCos(t) +p 


(6-37) 


I.  (6-35)  further  simplifies  to 


k N-1 

R(k+l)=-l-(-l)  P + 

l-pCos(t)+p  Cos(N-l)t-p  Cos(Nt) 

2- 2 

l-2pCos(t)+p 


(6-38) 


where  t=-~ k.  Equation  (6-38)  is  the  desired  result  which 
provides  a closed  form  solution  for  the  eigenvalues  of  the 
matrix  in  fig.  (6-4).  In  order  to  use  the  matrix  in 
fig.  (6-4)  in  the  fast  Weiner  filter  equation,  fig.  (6-4) 
must  represent  a covariance  matrix  for  real  data.  Thus, 
the  matrix  of  fig.  (6-4)  is  required  to  have  nonnegative 
eigenvalues,  which  is  the  subject  of  the  following  lemma. 


Lemma  (6-1).  Figure  (6-4)  illustrates  a nonnegative 
definite  matrix  which  becomes  singular  only  when  p=l. 
Hence,  this  matrix  is  positive  definite  for  all  the  values 

of  p such  that  0£p<l. 

Proof.  It  must  be  illustrated  here  that  R(k+1)  of 


the  increments  of  k. 


eq.  (6-38)  is  nonnegative  for  all 

n 

Since  t=-^;^k,  then 

(N-1) t=  n k 


(6-39) 


and 


Cos[(N-l)t]=(-l)^ 


(6-40) 


Also 


Nt=  n k+ 


n 

N-l 


k 


(6-41) 


thus 


Cos(Nt)  = (-l)^Cos(t)  (g_42) 

Equations  (6-40)  and  (6-42)  can  be  substituted  in 
eq.  (6-38)  to  give 


k N-l  l-pCos(t)  + (-l)^p^ 
R(k+1 ) =-l+ (-1 ) p +2 


, -1 

-(-1)  p Cos(t) 


l-2p+Cos  ( t)  +p“^ 


or 


(6-43) 


R(k+i;= 


(1-p^  [i-(-i)^p^-l] 
l-2pCos ( t) +p 


(6-44) 
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since  P<1,  the  numerator  of  the  above  equation  U always 
positive,  and  also  since  Cos(t)<l  the  denominator  is  also 
always  positive,  thus  R(k+1)  is  larger  than  zero  for  all 
values  of  k and  0<P<1.  If  P=l- 

the  exception  of  the  case  when  t is  the  zero  angle  (k=0) . 
If  t=0 , then  Cos(t)=l,  and 


R(l)  = 


M-1 

(1-p)  (l-*-p)(l-P  ) 

U-P? 


(6-45) 


or 


R(l)  = 


N-1 

(1+p) (1-P  ) 

1-P 


(6-46) 


NOW  if  P approaches  unity,  then  the  following  nonzero  valu. 

for  R(l)  is  obtained 

N-2 

OM  9 (6-47) 

R(l)  = l+P  p =2N-2 

pi 


Since  the  assumption  of  p<l  is  always  valid,  lemma  (6-1) 
has  illustrated  that  matrix  of  fig.  (6-4)  is  a nonsingular 

covariance  matrix. 

The  developments  of  the  fast  Wiener  estimator  are 
summarized  in  the  following  steps 
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1.  The  circular  image  degradation  model  is 
defined  by  equation  (6-16). 

2.  An  extended  object  vector  is  defined  by 
augmenting  a vector  ^ certain  statistical 
characteristics  to  the  actual  object  vector  f. 
Equation  (6-18)  illustrates  this  vector. 

3.  A circulant  nonsingular  covariance  matrix  of 
form  fig.  (6-4)  is  constructed  to  be  used  in  the 
Fast  Weiner  filter  of  eq.  (6-17) 

4.  The  augmented  vector  y is  generated  - this 
step  is  yet  to  be  established. 

5.  The  observation  is  modified  to  correspond 
to  the  extended  object  -this  step  is  yet  to  be 
established. 

Continuing  from  step  4 of  above,  the  mean  zero  randoiii 
process  vector  ^ has  the  following  statistics 


E(ZZ^)  =C 


(6-48) 


and 


E(fx  )=C 


fy 


(6-49) 


where  is  the  lower  right  (N-2)x(N-2)  portion  of  matrix 
fig.  (6-4)  and  is  the  upper  right  Nx(N-2)  portion  of 


103 


••'TT  » '■»)« 


the  same  matrix.  Thus,  C,^  is  a Markovian  covariance  matrix 
of  size  (N-2)x(N-2),  and  C,  is  given  as 


N-2„N-3 
p P 

N-1  N-2 
P P • • 

N-2  N-1 
P P • ' 


2 

P P 

3 2 
P P 

4 3 
P P 


C = 


(6-50) 


3 4 

P P 

2 3 

P P 

2 

P P 


N- 1 N-2 
P P 
n-2  n-1 
P P 

N-3  N-1 
? P 


The  statistical  information  on  ^ can  be  employed  to 
generate  this  vector.  To  proceed,  let  the  entries  of  y be 
equal  to  a linear  combination  of  the  entries  of  the 
observed  physical  image  samples,  g,  plus  a noise  term. 
Assuming  Q represents  the  linear  operation,  vector  y can  be 
expressed  as 


y=^+u 


(6-51) 


where  u is  an  independent  noise  term.  The  observation 
vector  g is  given  by  eq.  (3-41).  Thus, 


y=QBf-tQn-tu 


(6-52) 


;.n 


r 


T T T 
Q -QVQ 


(6-58) 


F- 


ft' 

•ft 


Equation  (6-58)  is  used  to  generate  the  mean-zero 
independent  noise  term  u.  The  numerical  value  of  u must  be 
Summed  with  vector  ^ to  generate  vector  y in  order  to 
complete  step  4.  In  the  last  step,  the  extended 
observation  vector  £ must  be  established.  To  proceed  at 
this  point,  two  hypothetical  K— dimensional  observation 
vectors  and  ^^are  defined  as  follows:  let  be  formed 
by  augmenting  £ with  a vector  of  zeros.  The  resulting 
vector  is 


2i 


(6-59) 


And  £2  is  formed  by  artificially  degrading  (blurring)  a 
vector  of  zeros  augmented  by  y-  Vector  g^  has  the 
following  form 


£ =C 
^2  - 


(6-50) 


The  vector  is  found  as  the  sum  of  £j  and  £2 
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Hence , 


or 


(6-61b) 


3^.  +n 


(6-61C) 


where  n is  a noise  term  which  incorporates  the  original 

observed  noise  plus  the  uncertainty  introduced  by  a 

possible  windowing  operation  applied  on  the  observation. 

Equation  (6-61c)  represents  a circular  degradation 

phenomenon,  whose  object  vector  f has  a circulant 

nonsingular  covariance  matrix.  The  minimum  mean-square 

error  estimate  of  the  object,  f , is  given  as 

— c 

A -1-1  T -1  T 

f =(S  C +C  C)  C q (6-62) 

-c  -fc ""c 


where  S is  the  signal  to  noise  ratio.  Extraction  of  the 

A 

first  N entries  of  f results  in  the  true  object  estimate 


f.  Thus 


(6-63) 


^ N 
l=iiKl 


where  Slj,  is  the  proper  selection  matrix, 


6.4  A Comment  on  the  Optimality  of  the  Fast  Filter 

The  previous  section  developed  a fast  image 
restoration  technique  which  is  optimal  in  a minimum 
mean-square  error  sense.  It  should  be  noted,  however,  that 
the  optimality  of  this  filter  depends  greatly  on  the 
validity  of  the  statistical  assumptions  imposed  upon  the 
hypothetical  object  vector  f ^ . Although  the  proper 
selection  of  the  physical  samples  of  the  auxiliary  random 
process  vector  y is  essential  for  recovery  of  the  object 
with  the  least  mean-square  error,  the  error  analysis  of  the 
next  section  illustrates  that  the  increase  in  the  error  can 
be  extremely  small. 

Because  of  the  circularity  of  the  covariance  matrix 
Cj.  , the  few  end  pixels  of  the  true  object  f happen  to  be 
highly  correlated  with  certain  entries  of  y.  To  illustrate 
this  claim,  fig.  (6-5a)  shows  the  pixels  of  f arranged 
around  a circle.  Under  a Markovian  assumption,  it  is  noted 
that 


if- 


Elf  Y 1=  P 
N+l-i  1 


i<  ^ 
2 


i>  -N- 


(6-64) 


Equation  (6  64)  illustrates  that  the  central  pixels  of  f, 
where  the  subscript  i is  not  close  to  either  one  cr  N,  are 
just  slightly  correlated  with  any  pixel  in  On  the  other 
hand,  the  pixels  of  f lying  close  to  the  two  ends,  where  i 
is  either  close  to  one  or  N,  are  strongly  correlated  with 
y or  y^ 


1 


N-2 


In  view  of  the  above  argument,  a poor  estimation  of 
the  physical  samples  of  i would  mostly  affect  the  few 
pixels  of  the  object  estimate  at  its  two  ends.  Therefore, 
it  can  be  expected  that  a reasonably  small  amount  of  error 
in  the  vector  y would  not  influence  the  center  portion  of 
the  object  estimate  f.  Thus  it  appears  that  computation  of 
the  physical  samples  of  the  random  process  ^ can  be  avoided 
if  this  vector  is  replaced  by  its  uean.  In  this  case, 
since  ^ is  mean  zero,  the  hypothetical  object  vector  f^  has 
the  following  form 


f c = 


f 

0 


(6-65a) 


and  the  observation  is  given  by 


no 


(6-65b) 


Representing  the  fast  filter  by  U,  the  object  estimate  f is 

~ -c 

obtained  as 


A 

4 = (6-66) 


The  first  N entries  of  f result  in  the  true  object 

A ^ 

estimate  f.  Thus, 


A 

f 


N A 
SI  f 
— K ~c 


N 

= S1  Uq 
■~K  — c 


(6-67) 


It  should  be  noted  that  the  above  equation  represents  a 
suboptimal  estimate  of  the  object  f,  but  the  estimation 
technique  associated  with  eq.  (6-67)  illustrates  extreme 
amount  of  computational  efficiency.  The  efficiency  of  this 
technique  is  due  to  the  fact  that  only  Fourier  transform 
operation  is  employed  for  obtaining  the  estimate  f of 
eq.  (6-67).  The  next  section  illustrates  that  replacing 
vector  ^ by  its  mean  results  only  in  an  slight  increase  of 
the  error  in  the  estimation  of  the  few  pixels  lying  at  the 
two  ends  of  the  object  vector. 


6.  5 Error  Analysis 

In  order  to  analyze  the  error  for  both  ovardetermined 
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and  underdetermined  systems  usinq  a uniform  terminology, 
the  following  approach  has  been  adopted.  Assume  x is  an 
object  of  size  0 and  that  the  N middle  pixels  of  this 
object,  f,  are  to  be  estimated.  Vector  x is  blurred  to 
generate  the  observation  g.  In  fact,  2 consists  of  the 

blurred  object  plus  an  additive  noise  term  n.  Next,  the 
observation  g undergoes  the  windowing  operation  W.  This 
modified  observation  is  referred  to  as  g.  Figure  (6-5^) 
illustrates  the  vectors  and  the  corresponding  operations. 
Both  g and  g are  of  size  M,  where 


M=Q-L+1 


(6-68a) 


and 


N=M-L+1 


(6-68b) 


The  observation  2 is  obtained  from  the  following  expression 
(see  Sec.  3.  5) 


g=Bx4-n 


(6-69) 


where  B is  an  M by  Q blur  matrix  (underdetermined),  and  n 
is  the  white  noise  term.  The  output  of  the  windowing 

operation  is  given  by 
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r 

q=WBx+n 

To  estimate  the  object  f,  by  means  of  the  fast  Wiener 
technique  and  by  use  of  the  observation  q,  the  random 
process  vector  y must  be  obtained.  For  the  sake  of 
computational  simplicity,  this  vector  can  be  replaced  by 
its  mean  value  which  is  zero.  It  must  be  realized  that 
this  approximation  and  application  of  the  windowing 
operator  are  the  only  approximations  made  in  computing  the 
classical  Wiener  filter.  In  other  words,  if  the  vector  y 
is  generated  by  executing  the  process  explained  in 

j 

Sec.  (6.  3),  the  error  term  for  the  fast  filter  must  be 
identical  to  the  error  function  derived  for  the  classical 
Wiener  filter  (assuming  a black  background  for  the  scene, 
or  absence  of  blur).  As  will  be  illustrated  later,  the 
approximation  of  y by  its  mean  results  in  a small  increase 
in  the  variance  of  the  error.  Nevertheless,  considering 
the  computational  savings  the  assumption  yields,  this 
approximation  is  worthwhile.  Thus  the  estimate  f is 
obtained  by 


(6-70) 


where  U represents  the  fast  filter.  For  the  sake  of 
shortening  the  length  of  the  preceeding  equations,  let 


N NT 

) (6-72) 

K i\ 

where  T is  an  N by  M matrix  (filter).  Using  eq.  (6-72)  in 

eq.  (6-71b)  results  in 

A 

f=la  (6-73) 

Next,  2 can  be  substituted  for  as  follows 
A 

f=TWBxtTWn  (6-74) 

The  error  covariance  matrix  is  defined  as 


C^=E(  (f-f)  (f-f)”^  (6-75a) 
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simplest  term  E(fjE^)  , and  assuming  that  the  object  is  a 
sample  from  a Markov  process,  this  term  can  be  represented 
by  an  N by  N Markovian  covariance  matrix.  The  term  E(^  ) 
is  obtained  as 

AT*  T 

E(ff  )=E[ (TWBx+TWn) f 1 (6-76a) 


or 


AT  T 

E(ff  ) =WB[E  ) 1 


(6-76b) 


where  E(xf"^)  is  the  cross  covariance  of  a Q dimensional 
object  with  its  own  N middle  p.xels.  Presenting  this  term 

by 


E(ff'^)*TWBC, 


£bcf 


(6-77) 


whei 


re  Cj^fis  the  Q by  N center  portion  of  a Q by  Q Markovian 


trix.  The  term  E[^)  is  simply  obtained  by 


covariance  ma 


transposing  E(^  ).  Thus, 


T T 

E(^  ) = (TWBC^) 


(6-78) 


AAT  . . . , 

And  the  term  E(«  ) is  obtained  as 


AAT 


E(ff  ) =E  (TWBx-t-TWn)  (TWBx-tTWn) 


(6-79) 
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aat  t t t t t 

E ( ff  ) =TWBC  B W T +TWVW  T 


(6-80) 


where  is  a Q by  Q Markovian  covariance  matrix  and  V is 
the  noise  covariance.  Finally,  the  error  covariance  is 
given  by  the  following  equation 


C =TWBC  B^W^T^+Cf +TWVW^T^- ( TWBC  +(TWBC  ^ )^1 
— e ^x—  — — -I  — xf XI 


(6-81) 


Figure  (6-6)  illustrates  the  error  variance  for  an  object 
estimate  containing  129  pixels.  The  error  is  plot*-ed  for 
two  different  blurring  effects.  Plot  (a)  represents  the 
error  for  Gaussian-shape  blur  of  standard  deviation  2. 
Plot  (b)  shows  the  error  for  motion  blur.  In  both  cases  L 
is  15,  the  signal-to-noise  ratio  is  10,  and  the  element 
correlation  coefficient  is  0.95.  Notice  that  the  error  for 
the  first  and  the  last  few  pixels  is  considerably  higher 
than  the  remainder.  This  can  be  explained  by  the 
approximation  made  on  vector  j.  The  circularity  property 
of  the  current  model  assumes  a strong  correlation  between  ^ 
and  the  first  and  the  last  few  pixels  of  the  object.  As  ! 
consequence  of  approximating  ^ by  a vector  of  zeros,  these 
pixels  become  correlated  with  the  wrong  data  — the  zeros -- 
which  results  in  a higher  error  variance.  Figure  (6^7) 
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illustrates  the  error  for  a higher  element  correlation  of 
0.99.  As  expected,  both  error  curves  decrease 
considerably.  Figure  (6-8)  contains  two  plots  for  a 
Gaussian-shape  blur  of  the  same  standard  deviation  but 
having  different  element  correlations.  Plot  (a)  is 
obtained  for  a correlation  coefficient  of  0.95,  and  (b)  is 
obtained  by  assuming  the  coefficient  is  0.99.  Figure  (6-9) 
is  the  counterpart  of  the  previous  figure  for  motion  blur. 


If  an  underdetermined  system  is  used  to  model  the 
degradation  process,  the  object  estimate  x can  be  obtained 
using  the  classical  Wiener  filter.  The  estimate  x is  given 


A 


(6-82) 


where  U represents  the  Wiener  filter.  The  error  term  is 
defined  as 


A 

e=x-x 


(6-83) 


and  the  error  vector  for  the  N middle  pixels  of  x is  given 


A 


(6-84) 


where  S2„  is  the  appropriate  selection  matrix. 


P 


covariance  matrix  of  e is  obtained  as  follows.  Let  C 

~e 

represent  this  covariance  matrix,  so  that 


C 

~e 


, A aT 
= E[ (x-x) (x-x)  ] 


(6-85) 


or 


^ T aaT  a T 

C^=E{>^  )+E{)^  )-E{^  )-E{3«  ) (6-86) 

T 

wnere  E{)^  ) is  a Q by  Q Markovian  matrix,  and  E{xx  ) is 
obtained  from 


aaT  T T 

E{)n(  )=U{E(2£  ) )u 


Recall  that 

T T 

E(22  ) = E(^+n)  (Bx+n) 


(6~88) 


or 


T T 

E{22  )=^B  +V  (g_89j 

Equation  (6-89)  can  then  be  used  in  equation  (6-87)  to 
yield  the  following  equation 


aaT  T T T 

E(xx  )=UBC^  U +UVU 


(6-90) 
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in  a similar  manner,  E(>^  ) can  be  derived  as 


aT  t t 
)=C^B  U 


(6-91) 


and 


A 


E(xx  )=UBC 


(6-92) 


Equations  (6-90),  (6-91),  and  (6-92)  can  be  used  in 

eq.  (6-86)  to  give  the  final  expression  for  the  error 

variance 


Cg’C^+UBC^B  U +UVU  -C^B  U 


(6-93) 


The  error  covariance  matrix  for  the  N middle  pixels  of  the 
object  can  be  obtained  by  observing  eq.  (6-84) 


A T N N. 

E(f-f)  (f-f)  =S2^(S2q) 


(6-94) 


Equation  (6-94)  is  the  error  equation  for  the  N-dimensional 
object  f.  Notice  that,  unlike  the  previous  case 
(overdetermined  system),  the  whole  object,  x»  is  estimated 

XT 


first  and  then,  using  a selection  matrix,  the 

N-dimensional  object  estimate,  f,  is  obtained.  Figure 


(6-10)  illustrates  the  error  for  the  case  when  N is  17 , L 
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Figure  (6-10)  Mean- square  error  comparison  of  Wiener  filter  and  fast  computational 
Wiener  filter. 


equals  9,  and  the  number  of  observed  pixels,  is  25. 

Plot  (a)  shows  the  error  for  the  classical  Wiener  filter, 
and  (b)  is  the  error  plot  for  the  fast  filter.  in  this 
figure,  except  for  the  first  and  the  last  few  pixels,  the 
error  is  the  same  for  both  the  Wiener  and  the  fast  Wiener 
filters.  As  was  pointed  out  earlier,  this  phenomenon  is 
caused  by  the  approximation  made  on  the  vector  in  both 
plots  the  blur  is  of  Gaussian  shape  with  a standard 
deviation  of  1,  and  the  correlation  coefficient  is  assumed 
to  be  0.7.  The  signal-to-noise  ratio  is  5.  Figure  (6-11) 
illustrates  the  error  when  there  is  no  blur  and  the  image 
degradation  is  due  only  to  additive  white  noise.  The 
signal-to-noise  ratio  is  assumed  to  be  10,  and  the  element 
correlation  to  be  zero.  In  this  example,  both  error 
functions  are  equal.  Since  no  correlation  between  the 
elements  is  assumed,  the  approximation  made  on  vector  y 
does  not  affect  the  error  plots.  Also,  in  the  absence  of 
any  blur  , the  windowing  operation  does  not  introduce  any 
uncertainty.  Figure  (6-12)  illustrates  the  error  curves 
for  a Gaussian-shape  blur  having  standard  deviation  of  2, 
and  an  element  correlation  of  0.95.  The  signal-to-noise 
ratio  is  assumed  to  be  100. 

6.  6 Experimental  Results 

To  illustrate  the  function  of  the  fast  filter. 
Fig.  (6-13a)  is  selected  as  a test  scene.  This  scene  is 
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f (a)  Original 


(c)  Restored 


(b)  Observed 


with  each 


represented  by  an  array  of  256  by  256  pixels, 
pixel  value  uniformly  quantized  to  256  levels.  This  image 
was  displayed  on  a flying  spot  scanner  cathode  ray  tube 
display  and  photographed  with  polaroid  type  52  film. 
Figure  (6-13b)  represents  the  image  after  it  undergoesing  a 
motion-blur  degradation  having  an  impulse  response  length 
of  15  pixels  (L=15).  To  observe  the  work  of  the  filter  on 
an  object  with  unknown  background,  only  the  middle  N 
( N <256  ) pixels  have  been  restored.  The  restored  object, 
in  this  example,  contains  129  (N=129)  pixels,  and  it  is 
placed  in  the  blurred  background  to  make  the  comparison 
between  the  observation  and  the  restored  image  simple,  as 
fig.  (6-130  shows,  in  this  figure,  the  background  of  the 
restored  object  is  obtained  by  carefully  extracting  the 
appropriate  section  of  the  blurry  observation.  This  kind 
of  object  and  background  combination  enables  the  observer 
to  easily  view  the  improvement  made  on  the  center  region  of 
the  scene.  Hence,  the  image  in  fig.  (6-13c)  can  be 

interpreted  as  the  blurry  observation  with  an  enhancing 
aperture  placed  in  front  of  the  scene.  In  this  particular 
example  the  hypothetical  aperture  contains  256  by  129 
pixels.  Since  motion  blur  is  a one  dimensional 
degradation,  each  line  of  the  observed  image  has  been 
processed  separately.  At  the  first  step,  143 

(143=129+15-1)  pixels  of  a line  of  the  observed  image  have 
been  extracted.  Next,  these  pixels,  after  undergoing  a 
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windowing  operation,  have  been  augmented  by  a vector  of 
zeros  to  form  a hypothetical  observation  vector  of  size  256 
(Note  that  the  size  of  the  hypothetical  observed  vector  g 
is  determined  by  eq.  (6-29) . After  the  hypothetical 
observation  is  formed,  it  is  Fourier  transformed,  and  then 
this  transformed  vector  is  multiplied  by  the  corresponding 
filter  coefficients  in  the  Fourier  domain  (a  scalar 

I operation) . Inverse  Fourier  transforming  the  filtered 

I vector  generates  the  hypothetical  object  vector  f . The 

first  129  (N=129)  pixels  of  this  vector  form  the  true 

object  vector 

6.  7 The  Problem  of  Unknown  Point  Spread  Function 

The  process  of  restoring  images  when  the  point  spread 
function  of  the  degrading  phenomenon  is  not  k lown  a priori 
is  usually  referred  to  as  blind  deconvolution  [6-10] . 
This  kind  of  problem  arises  when  the  characteristics  of  the 
imaging  system  are  not  known  to  the  observer,  and  thus  the 
impulse  response  must  be  directly  measured  from  the 
observed  image.  In  theory,  the  point  spread  function  can 
be  simply  obtained  by  a direct  measurement  of  the  image 
that  resuls  from  a point  source  of  light.  Such  an 

experimental  computation  of  the  point  spread  function  is 
severely  limited  in  practice  because  of  the  lack  of  real 
point  sources  in  the  original  scene.  A similar  technique 
known  as  edge  measurement  is  an  alternative  choice  when  the 
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point  spread  function  is  isotropic.  Unlike  isolated  point 
sources  of  light  which  rarely  occur  in  a scene,  edges  are 


abundant  in  most  images.  To  illustrate  how  an  edge 
measurement  can  help  to  determine  the  point  spread 
function,  assume  U(x)  represents  an  object  function  of  the 
following  form 

x>0 

x=0  (6-95) 
x<0 

The  function  U(x)  is  known  as  the  unit  step  function.  It 
is  a widely  accepted  discipline  to  represent  a point  object 
by  the  derivative  of  the  unit  step  function  [6-11],  [6-12]. 
Thus 

^(x)  = iV(x)  (6-96) 
dx 

where  «i(x)  is  now  the  mathematical  notation  for  a point 
object.  Assume  that  the  object  U(x)  has  undergone  a 
degradation  with  impulse  response  h(x).  Note  that  this 
function,  h(x),  is  yet  to  be  determined,  of  course.  The 
observation  g(x)  is  given  by  the  following  equation 


The  limits  of  the  above  definite  integral  are  omitted  since 
they  plsy  rule  in  this  analysis.  Differentiating  the 
above  equation  gives  rise  to  the  interesting  resultation 


or 


[-i-U(x-t)  ]h(t)dt 
dx 


A (x-t)h(t)dt 


(f-98) 


(6-99) 


notice  that  the  right  hand  side  of  the  above  equation  (by 
definition)  equals  the  impulse  response  h(x).  Hence 


dx 

Equation  (6-100)  implies  that  differentiating  the  image 
pattern  associated  with  an  edge  is  equivalent  to 
determination  of  the  impulse  response  of  a linear 
shift-invariant  optical  system.  Figure  (6-14a)  contains  an 
image  photographed  by  a SEM  Cambridge  stereo  scan  type 
S4-10  electron  microscope.  The  image  represents  a Ferrite 
(iron)  particle  taken  from  the  record  side  of  an  audio 
magnetic  tape.  The  magnification  ratio  is  130,000  to  1. 
Since  this  amount  of  magnification  equals  the  limiting 
power  of  the  system,  the  resulting  image  is  blurry.  To 
assure  a nondegraded  image,  it  has  been  experienced  that 
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the  magnification  ratio  on  this  system  must  stay  below 
100,000.  Fortunately  the  image  contains  many  edges  which 
can  help  to  determine  the  impulse  response  of  the  electron 


microscope.  The  estimated  impulse  response  in  this  example* 
was  approximated  by  a separable  two-dimensional  Gaussian 
function  of  standard  deviation  2.  Figure  (6-14b) 
illustrates  the  same  image  after  the  center  section  has 
been  restored.  The  size  of  the  restored  region  is  129x129, 
and  the  observed  image,  fig.  (6-14a) , contains  256x256 
pixels  which  have  been  uniformely  quantized  to  8 bits.  The 
restoration  technique  was  the  fast  Wiener  filtering  of 
Sec.  (6.  3).  Figure  (6-14c)  shows  the  crystal  after  its 
upper  left  corner  was  filtered  using  the  fast  filter. 
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7.  THE  PROBLEM  OF  POSITIVE  RESTORATION,  SUGGESTIONS 

For  further  research,  and  conclusions 

The  Wiener  restoration  approach  introduced  in  the 
previous  chapter  neglects  certain  a priori  information 
concerning  the  pictorial  data  which,  if  utilized  properly, 
could  improve  the  quality  of  the  restored  images. 

Non-negativeness  of  an  optical  scene,  for  instance,  is  a 
restriction  which  can  be  utilized  to  reduce  the  estimation 
error  variance.  Also,  it  has  been  shown  that  .he  visual 
quality  of  a restored  image  can  be  enhanced  if  the  human 
visual  sys'.tem  response  is  utilized  in  the  restoration 
process  [7-1].  The  following  section  summarizes  certain 

non-negative  restoration  approaches  which  are  adaptable  to 
the  fast  restoration  technique. 

7.  1 Constrained  Restoration 

Positive  restoration  refers  to  image  filtering 
techniques  which  employ  a positiveness  constraint  to 
improve  the  restoration  of  degraded  images.  It  has  been 
illustrated  that,  in  general,  linear  inequality  constraints 
reduce  the  error  covariance  of  the  object  estimate,  and 
also  improve  the  stablity  of  the  system  representing  the 
degradation  phenomenon.  Reference  [7-2]  illustrates  this 
claim  by  adopting  a numerical  analysis  approach  to  the 
problem,  although  the  increased  computational  requirements 
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limitation  on  the 


of  this  approach  impose  a severe 
allowable  sii;e  of  the  image. 

An  alternative  approach  to  this  problem  is  to  utilize 
the  Fourier  domain  properties  of  positive  signals.  Lukosz 
[7-3]  has  determined  upper  bounds  on  the  Fourier  pattern  of 
a non-negative  signal  and  has  shown  that  the  amplitude  of 
the  Fourier  transform  of  a non  negative,  band-limited 
signal  satisfies  the  constraint 


IG 


i lul  I i Ivl  I 
(u,v)  1<G(0,0)  1-  — ^ i 


(7-1) 


where  G(u,v)  represents  the  Fourier  transform  of  the 
non-negative  signal  g(x,y),  and  U and  V are  the  cutoff 
frequencies.  Figure  (7-la)  illustrates  this  bound.  It 
should  be  noted  that  eg.  (7-1)  is  a necessary  condition, 
but  not  necessarily  a sufficient  restriction.  7n  other 
words,  there  are  many  signals  which  satisfy  ineauality 
(7-1),  but  are  not  necessarily  non-negative.  The  Lukosz 

bound  in  its  present  form  does  not  apply  to  discrete 
signals,  and  the  only  constraint  which  appears  to  hold  for 
discrete  signals  can  be  expressed  as 


1 G ( i , j ) 1 <G ( 0 , 0 ) 


(7-2) 


where  the  G(i,j)  are  the  discrete  values  of  the  Fourier 
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transform  of  the  discrete  signal  q(l,m).  Figure  (7  lb) 
illustrates  this  bound. 

Constraints  having  the  form  shown  in  eq.  (7-1)  or 
eq.  (7-2)  are  very  simple  to  implement,  particularly,  if  in 
the  filtering  process,  the  Fourier  transforms  can  be 
computed  in  advance.  For  example,  implementation  of  the 
inequality  in  eq.  (7-2),  for  images  restored  by  the  fast 
Wiener  filter  of  Sec.  (6-3),  is  trivial  in  nature.  This  is 
because  the  Wiener  filtering  takes  place  in  the  Fourier 
domain,  and  all  that  remains  is  to  check  every  coefficient 
of  this  domain  against  the  d.c.  value  of  the  object, 

G(0,0).  If  any  coefficient  violates  this  bound,  the 
coefficient  will  be  automatically  decreased  until  the  bound 
is  satisfied. 

The  major  shortcoming  of  the  Fourier  domain  inequality 
restrictions  for  positive  image  restoration  stems  from  the 
fact  that  most  images  contain  small  quantities  of  high 
frequency  components,  and  thus  rarely  violate  an  inequality 
of  the  form  of  eq.  (7-2) . For  instance,  although  the 
restored  object  of  fig.  (6— 14b)  occasionally  contains 
negative  quantities,  it  never  violates  the  bound  described 
by  eq.  (7-2).  hence,  before  displaying  the  image,  the 
negative  entries  were  cliped  to  zero. 

A brute  force  approach  for  assuring  the  positiveness 
of  the  images  restored  by  the  fast  filter  can  be 


implemented  as  follows.  Considering  a one-dimensional 
case,  let  F(i)  -epresent  the  ith  Fourier  coefficient  of  a 
K-dimensional  image  estimate  f.  Starting  from  the 
coefficient  associated  with  the  lowest  freouency,  F(0),  an 
image  is  constructed  using  the  Fourier  basis  function 
corresponding  to  the  selected  coefficient,  its  complex 
conjugate,  and  all  the  lower  order  basis  vectors.  Note 
that  since  the  image  is  real,  every  Fourier  entry  is 
accompanied  by  another  coefficient  which  is  its  complex 
conjugate,  and  that  the  d.c.  term  F(0)  is  always  positive. 
Thus,  at  the  very  first  stage,  the  image  is  represented  by 
a gray  level  which  has  a numerical  value  equal  to  F(0);  at 
the  second  step,  the  two  next  lowest  frequency  vectors  are 
added.  This  process  is  continued  until  a negative  quantity 
is  detected  in  the  constructed  image  when,  in  this  case, 
the  two  coefficients  corresponding  to  the  very  last  two 
vectors  added  to  the  image  are  modified  to  retain  the  image 
positiveness.  Then,  the  procedure  is  carried  on  until  all 
the  coefficients  are  exh-iusted  and  the  image  is  totally 
created.  In  some  cases,  however,  the  reconstructed  image 
may  never  become  negative.  To  avoid  the  lengthy 
requirements  of  positive  restoration  in  such  a case,  the 
image  can  be  first  constructed  by  simply  inverse  Fourier 
transforming,  and  then  if  negative  quantities  are  found, 
the  positive  restoration  technique  described  above  can  be 
applied.  The  only  drawback  of  the  brute  force  method  is  in 
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its  computational  inefficiency.  It  takes  a considerable 
amount  of  computing  time  to  check  if  every  single  Fourier 
coefficient  retains  the  image  positiveness.  One  suggestion 
requiring  further  research  in  this  area  is  to  attempt  to 
establish  methods  to  predict  the  coefficients  which  are 
likely  to  produce  negative  quantities.  This  type  of 
coefficient  selection  is  likely  to  be  a function  of  the 
size  of  the  quantity  itself.  For  instance,  the  quantities 
which  are  very  small  need  not  be  checked  out,  for  they 
usually  cannot  give  rise  to  negative  entries  in  the  image 
vector;  even  if  they  produce  negative  quantities,  their 
modification  cannot  improve  the  pictorial  data 
substantially. 


7.  2 The  Fast  Wiener  Filter  and  the  Eye  Model 


It  is  known  that  a subjectively  optimal  image  estimate 

can  be  obtained  if  the  human  visual  system  characteristics 

are  employed  to  constrain  the  restoration  technique. 

Unfortunately,  because  of  the  inaccessibl ity  and  complexity 

of  the  visual  system,  the  true  nature  of  this  system  is  not 

fully  known.  However,  indirect  measurements  and  repetitive 

experiments  have  unveiled  some  of  the  mysteries  of  the 

process  of  vision.  For  instance,  for  low  contrast  images, 

the  frequency  response  of  the  visual  system  is  believed  to 

be  of  form  of  fig.  (7-2),  and,  it  is  known  that  the  human 

visual  system  responds  nonlinearly  to  incident  light 
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intensity  (7-4] . 


To  simplify  the  notation  let  S*  denote  an  operation 
modelling  the  human  visual  system.  For  instance  S*  can 
represent  a simple  logarithm  operation,  although,  in 
reality  S*  can  only  be  expressed  by  much  more  complicated 
functions.  The  complexity  of  the  eye  model  makes  it 
impractical  to  design  a minimum  mean-square  restoration 
technique  which  satisfies  the  constraint  defined  by  S*  . A 
good  approach  for  avoiding  this  hindrance  is  to  process  the 
image  signal  output  of  the  visual  system.  Let  g represent 
the  observed  signal,  and  let  q denote  the  same  observation 
after  going  through  the  system  S*.  It  should  be  noted  that 
q is  actually  the  image  which  is  perceived  by  the  brain. 
To  improve  q,  the  Wiener  filter  technique  can  be  applied  to 
minimize  the  observation  error  in  the  mean-square  sense. 
Figure  (7-3)  illustrates  this  approach.  Figure  (7-3a) 
represents  the  eye  model,  where  f denotes  the  object  signal 
and  d represents  the  signal  observed  by  the  brain.  Figure 
(7-3b)  illustrates  the  filtering  process.  In  this  figure 
S*  denotes  the  inverse  function  of  the  visual  system.  It 
should  'be  noted  that  S*  is  followed  immediately  by  the 
observers  eye  system  S*.  Thus,  the  final  result  is  shown 
in  fig.  (7-3c).  This  figure  illustrates  that,  by  Wiener 
filtering  the  signal  entering  the  brain  q,  the  complicated 
eye  model  constrait  can  be  avoided. 
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To  summarize  this  dissertation,  the  next  section 
briefly  reviews  this  work  and  suggests  means  of 
generalizing  some  of  the  techniques  introduced  in  the 
previous  chapters. 

7.  3 Extensions  to  the  fast  Wiener  filter 

It  was  shown  earlier  in  this  dissertation  that  the 
fast  Wiener  filter  in  the  frequency  domain  can  be  described 
by  the  following  scalar  equation 


a(n) 


W(n)  = 


a(n)  1^  +5"  %^(n) 


(7-3) 


where  a(n)  represents  the  n-th  eigenvalue  of  the 
corresponding  blur  matrix,  b(n)  represents  the  n-th 
eigenvalue  of  the  circulant  covariance  matrix,  and  S is  the 
signal-to-noise  ratio.  If  the  noise  power  approaches  zero, 
eq.  (7-3)  represents  the  inverse  filter 


I(n)  = 


(7-4) 


a(n) 


where  I(n)  represents  the  n-th  entry  of  the  inverse  filter 
in  the  frequency  domain. 


is 


It 

has  been 

shown 

that 

not 

usually 

the 

best 

the  minimum  mean-square  filter 
technique  for  image  restoration 
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applications.  This  is  merely  because  this  filter  does  not 
emphasize  high  frequency  components  of  an  image,  and  these 
components  are  very  important  for  visual  perception  [7-4]. 
Figure  (7-2)  shows  the  frequency  response  of  the  human 
visual  system,  illustrating  the  importance  of  high 
frequency  pictorial  information  for  a human  observer. 
Since  the  human  visual  system  attenuates  lower  , spatial 
frequencies,  the  higher  frequency  components  must  be  of 
greater  utility.  Considering  this  argument,  it  may  seem 
that  the  inverse  filter  approximates  the  human  visual 
system  closer  than  the  Wiener  filter.  Unfortunately,  in 
the  presence  of  noise,  the  inverse  filter  is  not  feasible. 
However,  a compromise  can  be  achieved  if  the  two  filters 
are  averaged  in  a geometrical  sense  [7-5].  A Geometrical 
Mean  filter  is  formulated  as 


s • 1 " s 

Va(n)/  \|  a(n)  j +S  b (n)/ 


(7-5) 


where  0<  s <1,  and  G(n)  is  the  n-th  entry  of  the  filter  in 
the  frequency  domain.  Clearly,  the  Wiener  and  the  inverse 
filters  ace  both  special  cases  of  the  Geometrical  Mean 
filter  for  s=0  and  s=l,  respectively.  Since  eq.  (7-5) 
defines  a scalar  operation,  this  filter  is  already  in  a 
computationally  efficient  (fast)  form.  This  holds  true 
since  b(n)  is  an  eigenvalue  of  the  circulant  covariance 

matrix.  To  apply  a fast  Geometrical  Mean  filter  to  a 
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blurry  observation,  the  same  steps  which  are  involved  in 
the  application  of  the  fast  Wiener  filter  must  be  observed. 
For  instance,  here  too  the  M-dimensional  observation  2. 
operated  upon  by  the  windowing  matrix  to  reduce  the 
wrap-around  error  and  then  is  placed  in  a longer  vector  of 
zeros  of  size  K,  where  K is  given  by  eq.  (6-29).  And, 
after  the  K-dimensional  observation  is  processed  by  the 
fast  flter,  a K-dimensional  object  estimate  results  in 
which  the  first  N entries  of  this  vector  are  the  true 
object  estimate  f.  Since  the  Geometrical  Mean  filter 
represents  a class  of  restoration  techniques  to  which  the 
Wiener  and  the  Inverse  filters  are  special  cases,  only  one 
general  software  (filter)  is  necessary  to  be  designed  to 
give  a wide  selection  of  restoration  techniques. 


7.  4 Summary  and  Conclusons 


This  dissertation  has 
efficient  image  restoration 
examples  have  been  presented  to 
of  these  techniques. 


developed 
techniques , 
illustrate 


computationally 
and  pictorial 
the  efficiency 


The  continuous  convolution  integral  has  been  utilized 
to  represent  linear  shift-invariant  degradation  phenomena. 
It  has  been  shown  that  the  discrete  image  degradation 
problem  can  be  modelled  by  an  over  determined  system,  if  the 
object  background  is  not  known.  To  avoid  this,  dual  model. 
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a simple  operation  has  been  introduced  which,  by  predicting 
the  object  background,  almost  eliminates  the  contribution 
of  the  background  to  the  pictorial  data.  Hence,  only  the 
overdetermined  model  is  employed  to  describe  the 
degradation  problem.  This  operation  is,  in  essence, 
designed  to  control  the  modelling  error.  In  the  absence  of 
noise,  the  inverse  filter  for  the  continuous  case  or  the 
matrix  pseudoinverse  for  the  discrete  data  can  be  emp?  lyed 
to  restore  degraded  images.  The  computational  shortcomings 
of  the  pseudoinverse  techniques  can  be  overcome  oy 
utilizing  the  Fourier  domain  properties  of  circulant 
matrices.  Simulated  pictorial  examples  were  used  to 
illustrate  this  point. 


In  a noisy  environment,  the  statistics  of  the  image 
and  the  noise  can  be  employed  to  control  high  frequency 
noise  oscillation  in  the  restored  data.  Based  on  this 
fact,  a minimum  mean-square  error  filter  can  be  constructed 
for  image  restoration.  Since  a filter  of  this  kind  is 
computationally  unattractive,  a fast  Wiener  filter  was 
introduced  for  noisy  image  restoration.  This  filter  was 
obtained  by  imposing  certain  modifications  on  the  observed 
image.  It  was  shown  that  the  observed  image  can  be 
operated  upon  to  modify  its  statistical  characteristics. 
Thus,  certain  operations  were  introduced  which,  when 
applied  rn  pictorial  data,  allow  the  modified  data  to  be 
characterized  statistically  by  a circulant  covariance 
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matrix.  Later,  these  operations  were  approximated  to  gain 
more  computational  speed  at  the  cost  of  a slight  increase 
in  the  mean-square  error.  An  error  study  then  illustrated 
that  for  most  of  the  length  of  the  object  vector,  the  error 
variance  associated  to  the  fast  Wiener  filter  is  equal  to 
the  one  associated  with  the  Wiener  filter  itself.  The 
slight  error  increase  occurs  only  at  the  beginning  and  the 
end  of  the  object  vector,  and  almost  vanishes  when  the 
element  correlation  becomes  small.  Unlike  the  classical 
Wiener  filter,  it  was  shown  that  the  fast  filter  is  capable 
of  operating  upon  large  images  without  the  need  to  break 
down  the  observation  into  small  blocks. 

Statistically  speaking,  it  was  assumed  that  a scene  is 
a sample  of  a Markovian  random  process.  This,  of  course, 
is  not  an  essential  restriction,  and  can  be  removed.  It 
appears  that  the  fast  filter  can  be  constructed  for  any 
image  which  possesses  a monotonically  decreasing 
correlation  function.  In  fact,  any  covariance  matrix  which 
can  be  extended  into  a larger,  circulant,  and  positive 
definite  matrix  characterizes  a random  image  that  can  be 
operated  by  the  fast  filter.  Thus,  lemma  (6-1)  can  be 
n:oved  for  certain  non-Mar kovian  sources  as  well. 

Since  the  main  objective  in  this  dissertation  has  been 
to  present  computationally  efficient  restoration 
techniques,  it  seems  proper  at  this  stage  to  point  out  that 
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the  classical  Wiener  approach  is  not  the  only  restoration 
technique  which  can  be  modified  for  computational 
efficiency.  In  fact,  as  previously  shown,  there  is  a class 
of  filters  which  can  be  efficiently  computed  through  a 
technique  similar  to  the  Wiener  approach.  Notice  that  the 
phrase  "computational  efficiency"  corresponds  to  both  the 
time  and  storage  requirements  of  a certain  restoration 
technique.  Hence,  a fast  counterpart  of  a restoration 
technique  is  computationally  advantageous,  in  the  sense 
that  the  fast  filter  permits  processing  of  large  images  in 
if0latively  small  amount  of  time. 


references 


1.  B.  Hunt,  "Digital  Image  Processing,"  g£.0£^  April 

1975,  pp.  693-708. 

2.  N.  Mascarenhas,  Digital  Image  Restoration  Unde£__a 
Regression  Model-  The  Unconstrained ^ing.a£ 

Inequality  Constrained  Approaches,  Ph . D.  Dissertation, 
University  of  Southern  California. 


3.  W.  Lukosz,  "Ubertragung  Micht-Negativer  Signale  durch 
Lineare  Filter,"  Optica  Octa,  9,  1962,  pp.  335-364. 


